Two factor production with triangle density that is greater zero everywhere: numerical example 1

c[q_,t_]=a*t*q+q^2/t-g*t;
cq[q_,t_]=a t+2q/t;
ct[q_,t_]=a q-q^2/t^2-g;
cqt[q_,t_]=a-2q/t^2;
cqqt[q_,t_]=-2/t^2;
cqtt[q_,t_]=4q/t^3;
tl=1/4; (*theta lower bar*)
tu=3/4;   (*theta upper bar*)
delta=tu-tl;
(*Assign some parameter values for a,b and g*)
b=16/10;(*principal's utility is bq *)
a=1;
g=1/3;
s[t_]=a t^2/2;

f[t_]=(2*9/10*(t-tl))/(delta^2)+1/(10*delta); (*density of theta*)
F[t_]=9/10(t-tl)^2/(delta^2)+(t-tl)/(10*delta);(*distribution*)

(*relaxed decision and optimal decision as function of eta*)
Assuming[tu≥t≥tl,Reduce[(b -cq[q,t])f[t]+(1-F[t]) cqt[q,t]==0,q,Reals]]
Assuming[tu≥t≥tl&&n≥0,Reduce[(b -cq[q,t])f[t]+(1-F[t]-n) cqt[q,t]==0,q,Reals]]

(*Straightforward calculation gives the first best and relaxed decision:*)

qfb[t_]=-(1/2) t (-b+a t);
qr[t_]=((-347 (t^2)+2944 (t^3)-2160 (t^4))/(330+1440 (t^2)));
pirt[t_]=-ct[qr[t],t];
pir[t_]= (∫_tl^t)pirt[s]\[DifferentialD]s;
Phir[t_,h_]=pir[t]-pir[h]-c[qr[h],h]+c[qr[h],t];

Plot[{qr[t],s[t],qfb[t]},{t,tl,tu}]

(*Here I check the sufficient condition in proposition 1. To do so, the mirror images q^v and q^s have to be computed. The results are stored as functions qs1 and qv1. The following Reduce command below verifies that qs1≤qv1 as the opposite statement is "False" in the relevant range. The following plot shows that the first best and the relaxed decision are above s(theta) and both are strictly monotonically increasing. Hence, assumption 1 is also satisfied.*)
VV[q_,t_]=(b q-c[q,t])f[t]+(1-F[t])ct[q,t];
Reduce[{ct[q,t]==ct[qs,t],qs>s[t],tl≤t≤tu,0≤q<s[t]},qs,Reals]
Reduce[{VV[qv,t]==VV[q,t],qv>s[t],0≤q<s[t],tl≤t≤tu},qv,Reals]

qs1[q_,t_]=t^2-q;
qv1[q_,t_]=((-165 q-347 (t^2)-720 q (t^2)+2944 (t^3)-2160 (t^4))/(165+720 (t^2)));
Reduce[{qs1[q,t]>qv1[q,t],0≤q<s[t],tl≤t≤tu},{q,t},Reals]

(*The following plot illustrates for which types the non-local incentive constraint is violated
under the relaxed solution*)

Plot3D[{Phir[t,h],0},{t,tl,tu},{h,tl,tu}]

Now take the parameter values above and use the algorithm. First step is to get for a given eta (i.e. "n") the corresponding decisions, profits and incentive constraints.

qn[t_,n_]=((-347 (t^2)-200 n (t^2)+2944 (t^3)-2160 (t^4))/(330-400 n+1440 (t^2)));
pin[t_,n_]=Assuming[{tl≤t≤tu,0≤n≤1},(∫_tl^t)-(qn[s,n]-(qn[s,n])^2/(s^2)-1/3)\[DifferentialD]s];
Phin[t_,h_,n_]=pin[t,n]-pin[h,n]-(h*qn[h,n]+qn[h,n]^2/h-1/3 h)+(t*qn[h,n]+qn[h,n]^2/t-1/3 t);

NMinimize[{Phin[t,h,n]/.{n->0},1/4≤h≤t≤3/4},{{t,7/10,3/4},{h,1/4,3/10}}]

Next, Phin is minimized and the minimizers as well as the minimum value is reported :

For[i=1,i<11,i++,Print[N[5*i/100,6]];Print[NMinimize[{Phin[t,h,n]/.{n->5*i/100},1/4≤h≤t≤3/4},{{t,7/10,3/4},{h,1/4,3/10}}]]]

(*The relevant \bar{eta} turns out to be around 0.47:*)

For[i=472980,i<472985,i++,Print[N[i/1000000,6]];Print[NMinimize[{Phin[t,h,n]/.{n->i/1000000},1/4≤h≤t≤3/4},{{t,7/10,3/4},{h,1/4,3/10}}]]]

(*Ok, eta bar would be 0.472982. The constraint is only binding from top to bottom and not at any interior types.*)

(*Now, the solution denoted by qs. Obviously, all t have the same eta, i.e. eta=0.472982*

qs[t_]=qn[t,47298/100000]

Plot[{qs[t],s[t],qfb[t],qr[t]},{t,tl,tu},PlotStyle->{Red,Gray,Green,Blue}]

(*Plotting the rent as functon of type*)

Plot[pin[t,472982/1000000],{t,1/4,3/4}]

FullSimplify[{qs[t],s[t],qfb[t],qr[t],qn[t,n]}]