%u(c)=c^s/s. p=probability of high state, h=endowment in high state,
%l=endowment in low state, d=delta.
% initialize needed variables and create initial values for coefficient
% matrices
global s p d h l xbar
s=-2;
p=.5;
d=.95;
h=1.9;
l=0.1;
%in the first pass, we use the complete information value function val1 as
%next period's value function; we create an approximation of this period's
%value function under incomplete information and iterate on that in the
%following passes. The functional form of val1 depends on whether d is large
%enough so that efficient risk sharing is possible under complete info.

klow=(p*(h^s)+(1-p)*l^s)/(s*(1-d));
options = optimset('LargeScale', 'off');
n=1;

if 1-(1/((d*p*(1+h)^s+(1-d*p)*(1+l)^s)^(1/s)))>=((((1-d*(1-p))*h^s+d*(1-p)*l^s)/((1-d*(1-p))*(1+h)^s+d*(1-p)*(1+l)^s)))^(1/s)
    shar=1;
    fid=fopen('val1.m', 'w');
    fprintf(fid, 'function F=val1(k)\n');
    fprintf(fid, '%expected value function under complete information \n');
    fprintf(fid, 'global s p d h l\n');
    fprintf(fid, 'if k<=((p*(1+h)^s+(1-p)*(1+l)^s)/(s*(1-d)))*(((1-d*(1-p))*h^s+d*(1-p)*l^s)/((1-d*(1-p))*(1+h)^s+d*(1-p)*(1+l)^s))\n');
    fprintf(fid, 'f=((p*(((1-d*(1-p))*(1+h)^s+d*(1-p)*(1+l)^s)^(1/s)-((1-d*(1-p))*h^s+d*(1-p)*l^s)^(1/s))^s+(((1-d)*(1-p))^(1/s)*(1+l)-(s*(1-d)*(1-d*(1-p))*k-p*((1-d*(1-p))*h^s+d*(1-p)*l^s))^(1/s))^s)/(s*(1-d*(1-p))*(1-d)));\n');
    fprintf(fid, 'fprime=(((((1-d)*(1-p))^(1/s)*(1+l))/((s*(1-d)*(1-d*(1-p))*k-p*((1-d*(1-p))*h^s+d*(1-p)*l^s))^(1/s)))-1)^(s-1);\n');
    fprintf(fid, 'elseif k<=((p*(1+h)^s+(1-p)*(1+l)^s)/(s*(1-d)))*((((d*p*(1+h)^s+(1-d*p)*(1+l)^s)^(1/s)-1)^s)/(d*p*(1+h)^s+(1-d*p)*(1+l)^s))\n');
    fprintf(fid, 'f=((((p*(1+h)^s+(1-p)*(1+l)^s)^(1/s)-((1-d)*s*k)^(1/s))^s)/(s*(1-d)));\n');
    fprintf(fid, 'fprime=((((p*(1+h)^s+(1-p)*(1+l)^s)^(1/s))/(((1-d)*s*k)^(1/s)))-1)^(s-1);\n');
    fprintf(fid, 'else\n');
    fprintf(fid, 'f=((((1+h)*((1-d)*p)^(1/s)-(s*(1-d)*(1-d*p)*k-(1-p)*((d*p*(1+h)^s+(1-d*p)*(1+l)^s)^(1/s)-1)^s)^(1/s))^s+(1-p))/(s*(1-d*p)*(1-d)));\n');
    fprintf(fid, 'fprime=(((((1-d)*p)^(1/s)*(1+h))/((s*(1-d)*(1-d*p)*k-(1-p)*((d*p*(1+h)^s+(1-d*p)*(1+l)^s)^(1/s)-1)^s)^(1/s)))-1)^(s-1);\n');
    fprintf(fid, 'end\n');
    fprintf(fid, 'F=[f;fprime];\n');
    fclose(fid);
    fid=fopen('vh1.m', 'w');
    fprintf(fid, 'function v=vh1(gam)\n');
    fprintf(fid, '%Value function in high state \n');
    fprintf(fid, 'global s p d h l\n');
    fprintf(fid, 'if gam<=(((1-d*(1-p))*(1+h)^s+d*(1-p)*(1+l)^s)/(s*(1-d)))*((((d*p*(1+h)^s+(1-d*p)*(1+l)^s)^(1/s)-1)^s)/(d*p*(1+h)^s+(1-d*p)*(1+l)^s))\n');
    fprintf(fid, 'Gx=(gam*(p*(1+h)^s+(1-p)*(1+l)^s))/((1-d*(1-p))*(1+h)^s+d*(1-p)*(1+l)^s);\n');
    fprintf(fid, 'gx=(s*gam-d*s*Gx)^(1/s);\n');
    fprintf(fid, 'else\n');
    fprintf(fid, 'Gx=((gam*(1-d)*p*s+(1-p)*((d*p*(1+h)^s+(1-d*p)*(1+l)^s)^(1/s)-1)^s)/(s*(1-d)));\n');
    fprintf(fid, 'gx=(s*gam-d*s*Gx)^(1/s);\n');
    fprintf(fid, 'end\n');
    fprintf(fid, 'fval=(1/s)*(1+h-gx)^s+d*[1 0]*val1(Gx);\n');
    fprintf(fid, 'vprime=((1+h-gx)/gx)^(s-1);\n');
    fprintf(fid, 'v=[gx; fval; vprime];\n');
    fclose(fid);
    khigh0=((p*(h^s)*(1-d)+(1-p)*((d*p*(1+h)^s+(1-d*p)*(1+l)^s)^(1/s)-1)^s)/(s*(1-d*p)*(1-d)));
    bhigh=(1+l)*(1-(1/((d*p*(1+h)^s+(1-d*p)*(1+l)^s)^(1/s))));
    Bhigh=((p*(1+h)^s+(1-p)*(1+l)^s)/(s*(1-d)))*((((d*p*(1+h)^s+(1-d*p)*(1+l)^s)^(1/s)-1)^s)/(d*p*(1+h)^s+(1-d*p)*(1+l)^s));
else
    shar=0;
    fid=fopen('autark.m', 'w');
    fprintf(fid, 'function [F,J]=autark(x)\n');
    fprintf(fid, '%system of equations defining compinfo solution,x=g,b \n');
    fprintf(fid, 'global s p d h l\n');
    fprintf(fid, 'F(1)=(1-d*(1-p))*x(1)^s+d*(1-p)*x(2)^s-(1-d*(1-p))*h^s-d*(1-p)*l^s;\n');
    fprintf(fid, 'F(2)=d*p*(1+h-x(1))^s+(1-d*p)*(1+l-x(2))^s-1;\n');
    fprintf(fid, 'if nargout>1\n');
    fprintf(fid, 'J=[s*(1-d*(1-p))*x(1)^(s-1) s*d*(1-p)*x(2)^(s-1); -s*d*p*(1+h-x(1))^(s-1) -s*(1-d*p)*(1+l-x(2))^(s-1)];\n');
    fprintf(fid, 'end\n');
    fclose(fid);
    options = optimset('Jacobian','on', 'LargeScale', 'off');
    x0=[1.2 .6];
    xbar=fsolve('autark', x0, options);
    fid=fopen('val1.m', 'w');
    fprintf(fid, 'function F=val1(k)\n');
    fprintf(fid, '%expected value function under complete information \n');
    fprintf(fid, 'global s p d h l xbar\n');
    fprintf(fid, 'if k<=((p*xbar(1)^s+(1-p)*xbar(2)^s)/(s*(1-d)))\n');
    fprintf(fid, 'f=((p*((1-d*(1-p))*(1+h-xbar(1))^s+d*(1-p)*(1+l-xbar(2)))+((1+l)*((1-p)*(1-d))^(1/s)-((1-d*(1-p))*(1-d)*s*k-p*((1-d*(1-p))*h^s+d*(1-p)*l^s))^(1/s))^s)/(s*(1-d*(1-p))*(1-d)));\n');
    fprintf(fid, 'fprime=(((((1-p)*(1-d))^(1/s)*(1+l))/((s*(1-d*(1-p))*(1-d)*k-p*((1-d*(1-p))*h^s+d*(1-p)*l^s))^(1/s)))-1)^(s-1);\n');
    fprintf(fid, 'else\n');
    fprintf(fid, 'f=((((1+h)*(p*(1-d))^(1/s)-((1-d*p)*(1-d)*s*k-(1-p)*(d*p*xbar(1)^s+(1-d*p)*xbar(2)^s))^(1/s))^s+(1-p))/(s*(1-d*p)*(1-d)));\n');
    fprintf(fid, 'fprime=((((p*(1-d))^(1/s)*(1+h))/((s*(1-d*p)*(1-d)*k-(1-p)*(d*p*xbar(1)^s+(1-d*p)*xbar(2)^s))^(1/s)))-1)^(s-1);\n');
    fprintf(fid, 'end\n');
    fprintf(fid, 'F=[f;fprime];\n');
    fclose(fid);
    fid=fopen('vh1.m', 'w');
    fprintf(fid, 'function v=vh1(gam)\n');
    fprintf(fid, '%Value function in high state \n');
    fprintf(fid, 'global s p d h l xbar\n');
    fprintf(fid, 'if gam>=(((1-d*(1-p))*xbar(1)^s+d*(1-p)*xbar(2)^s)/((1-d)*s))\n');
    fprintf(fid, 'Gx=((gam*p*(1-d)*s+(1-p)*(d*p*xbar(1)^s+(1-d*p)*xbar(2)^s))/(s*(1-d)));\n');
    fprintf(fid, 'gx=(s*gam-d*s*Gx)^(1/s);\n');
    fprintf(fid, 'else\n');
    fprintf(fid, 'Gx=((gam*(1-p)*(1-d)*s*(1+l)^s+p*((1-d*(1-p))*h^s+d*(1-p)*l^s)*(1+h)^s)/(s*(1-d*(1-p))*(1-d)*(1+h)^s+(1-p)*(1-d)*d*s*(1+l)^s));\n');
    fprintf(fid, 'gx=(s*gam-d*s*Gx)^(1/s);\n');
    fprintf(fid, 'end\n');
    fprintf(fid, 'fval=(1/s)*(1+h-gx)^s+d*[1 0]*val1(Gx);\n');
    fprintf(fid, 'vprime=((1+h-gx)/gx)^(s-1);\n');
    fprintf(fid, 'v=[gx; fval; vprime];\n');
    fclose(fid);
    options = optimset('Jacobian','off', 'LargeScale', 'off');
    khigh0=((p*(1-d)*h^s+(1-p)*(d*p*xbar(1)^s+(1-d*p)*xbar(2)^s))/((1-d*p)*(1-d)*s));
    fid=fopen('autfoc.m', 'w');
    fprintf(fid, 'function f=autfoc(B)\n');
    fprintf(fid, 'global s p d h l xbar\n');
    fprintf(fid, 'f=-(1/s)*(1+l-((1/(1-d))-d*s*[1 0]*val1(B))^(1/s))^s-d*B;\n');
    fclose(fid);
    Bhigh=fminbnd('autfoc', klow, khigh0, options);
    bhigh=1+l-((1/(1-d))-d*s*[1 0]*val1(Bhigh))^(1/s);
end
% to test that the value function does not dip below the goal, we need
% to know that the revelation constraint holds at the new khigh

if (h^s)/s - ((bhigh + h-l)^s)/s +d*(khigh0 - Bhigh) < 0
	error ('RH violated at khigh');
end

kr = (1-p)*(bhigh^s)/s + p*((bhigh+h-l)^s)/s + d*Bhigh;
khigh1 = p*((h^s)/s  + d*khigh0) +(1-p)*((bhigh^s)/s + d*Bhigh);
betahigh = (bhigh^s)/s + d*Bhigh;

% we now evaluate the new value function at 20 points; note that to find
% the proper curvature, points are more closely spaced at the rim.
% Solutions are given for khigh, klow and kr.

t(1) = klow;
hsol=vh1(((1-d*(1-p))*(h^s)+d*(1-p)*l^s)/(s*(1-d)));
w1(1) = p*[0 1 0]*hsol+(1-p)*(1/s + d*[1 0]*val1(klow));
bsol(1)=l;
q=1/(1+((1-p)/p)*(l/h)^(s-1));
di(1) = -q*[0 0 1]*hsol-(1-q)*(1/l)^(s-1);

t(22) = khigh1;
w1(23) = 1/(s*(1-d));
di(23) = -1/h^(s-1);
bsol(23) = bhigh;

for i = 2:11
	t(i) = klow + ((khigh1 - klow)/200)*(i - 1)^2;
end
for i=12:21
    t(i) = khigh1 - ((khigh1 - klow)/242)*(22-i)^2;
end
%specifying the solution for k>=kr
nr = min(find(t>= kr));
t(nr+1:23)=t(nr:22);
t(nr)=kr;
for i = nr:22
	hsol = vh1((t(i) -(1-p)*betahigh)/p);
	w1(i) = p*[0 1 0]*hsol + (1-p)/(s*(1-d));
	di(i) = - [0 0 1]*hsol;
    bsol(i) = bhigh;
end
% next we check whether P2 binds at some k
%for this, we need the function foctest1, which is the first order condition to W(k).
fid=fopen('foctest1.m', 'w');
fprintf(fid, 'function f = foctest1(b, k);\n');
fprintf(fid, 'global s p d h l\n');
fprintf(fid, 'f = p*[0 0 1]*vh1((k+(1-p)*((b+h-l)^s-b^s)/s))*((b^(s-1)-(b+h-l)^(s-1))/(p*(b+h-l)^(s-1)+(1-p)*b^(s-1)))-(1+l-b)^(s-1)*(1/(p*(b+h-l)^(s-1)+(1-p)*b^(s-1)))+[0 1]*val1((k-(p*(b+h-l)^s+(1-p)*b^s)/s)/d);\n');
fclose(fid);
%also, p2test1 tells us about 2's participation constraint: for any b and k
%at which the revelation constraint holds, p2test1 shows whether P2L is
%slack, tight or violated.
fid=fopen('p2test1.m', 'w');
fprintf(fid, 'function f = p2test1(b, k);\n');
fprintf(fid, 'global s p d h l xbar\n');
fprintf(fid, 'f = (1/s)*(1+l-b)^s+d*[1 0]*val1((k-(p*(b+h-l)^s+(1-p)*b^s)/s)/d)-1/(s*(1-d));\n');
fclose(fid);

for i = nr-1:(-1):2
    if p2test1(1+l-eps,t(i))>=0
        np=i;
        break
    else
        k=t(i);
        bsol(i) = fzero('p2test1', [bhigh 1+l-eps], options, k);
    end
    if foctest1(bsol(i),t(i)) > 0
		hsol = vh1(t(i)+(1-p)*((bsol(i)+h-l)^s-bsol(i)^s)/s);
        w1(i) = p*[0 1 0]*hsol+(1-p)/(s*(1-d));
        vB=[0 1]*val1((t(i)-(p*(bsol(i)+h-l)^s+(1-p)*bsol(i)^s)/s)/d);
		di(i) = -p*[0 0 1]*hsol*(1+((1-p)*(bsol(i)^(s-1)-(bsol(i)+h-l)^(s-1)))/((1+l-bsol(i))^(s-1)/vB-p*(bsol(i)+h-l)^(s-1)-(1-p)*bsol(i)^(s-1)));
	else
		np = i;
		break
	end
end
%We proceed using foctest to find the solution.
%Having this solution we check that indeed B>=klow; if not, we set B=klow.
%For this we need a function that (when set to zero) will give the correct b(k)
%if B(k)=klow, for any k.
fid=fopen('Blow.m', 'w');
fprintf(fid, 'function f = Blow(b, k,klow);\n');
fprintf(fid, 'global s p d h l\n');
fprintf(fid, 'f =k-d*klow-(p*(b+h-l)^s+(1-p)*b^s)/s;\n');
fclose(fid);
%also one that (when set to zero) gives the correct b if B=khigh, for any
%k.
fid=fopen('Bup.m', 'w');
fprintf(fid, 'function f = Bup(b, k, khigh0);\n');
fprintf(fid, 'global s p d h l\n');
fprintf(fid, 'f =k-d*khigh0-(p*(b+h-l)^s+(1-p)*b^s)/s;\n');
fclose(fid);

for i=2:np
    k=t(i);
    mick=0;
    if t(i)>= (p*(1+h-eps)^s+(1-p)*(1+l-eps)^s)/s+d*klow
        bone=1+l-eps;
    else
        bone=fzero('Blow', [l 1+l-eps], options, k,klow);
        mick=1;
    end
    if t(i)<=(p*h^s+(1-p)*l^s)/s+d*khigh0
        btwo=l;
    else
        btwo=fzero('Bup', [l 1+l-eps], options, k, khigh0);
    end
    if foctest1(bone,k)>=0
        bsol(i)=bone;
    elseif foctest1(btwo,k)<=0
        bsol(i)=btwo;
        mick=0;
    else
        bsol(i)=fzero('foctest1', [bone, btwo], options, k);
        mick=0;
    end
    hsol=vh1(t(i)+(1-p)*((bsol(i)+h-l)^s-bsol(i)^s)/s);
    if mick==1
        w1(i) = p*[0 1 0]*hsol+(1-p)*((1/s)*(1+l-bsol(i))^s+d*[1 0]*val1(klow));
    else
        w1(i) = p*[0 1 0]*hsol+(1-p)*((1/s)*(1+l-bsol(i))^s+d*[1 0]*val1((t(i)-(p*(bsol(i)+h-l)^s+(1-p)*bsol(i)^s)/s)/d));
    end
    q=1/(1+((1-p)/p)*(bsol(i)/(bsol(i)+h-l))^(s-1));
    di(i) = -q*[0 0 1]*hsol-(1-q)*((1+l-bsol(i))/bsol(i))^(s-1);
end

clear functions
    
% we now compute the spline
for i=1:23
    taug(2*i-1)=t(i);
end
for i = 2:23
	if di(i - 1) < di(i)
		error ('value function not concave');
	end
end
for i = 2:23
    z(i) = (w1(i)-w1(i-1))/(t(i)-t(i-1));
    if z(i-1)<z(i)
        error('value function not concave at all');
    end
end
	
for i = 2:23	
	if .5*(di(i) + di(i-1)) == z(i)
		taug(2*(i-1)) = (t(i) + t(i-1))/2;
		coeff(2*i-3,1)=(di(i)-di(i-1))/(2*(t(i)-t(i-1)));
        coeff(2*i-3,2)= di(i-1)-2*t(i-1)*coeff(2*i-3,1);
        coeff(2*i-3,3) =w1(i-1)-di(i-1)*t(i-1)+t(i-1)^2;
		coeff(2*i-2,1)=coeff(2*i-3,1);
        coeff(2*i-2,2)=coeff(2*i-3,2);
        coeff(2*i-2,3)=coeff(2*i-3,3);
	else
		if (di(i-1)-z(i))*(di(i)-z(i))>=0
            error ('going weird again');
			taug(2*(i-1))=(t(i)+t(i-1))/2;
		elseif abs(di(i-1)-z(i))> abs(di(i)-z(i))
			taug(2*(i-1))=t(i-1)+(t(i)-t(i-1))*(di(i)-z(i))/(di(i)-di(i-1));
		else
			taug(2*(i-1))=t(i)+(t(i)-t(i-1))*(di(i-1)-z(i))/(di(i)-di(i-1));
		end
        coeff(2*i-3,1) = (2*(w1(i) - w1(i-1))-di(i-1)*(taug(2*(i-1))-2*t(i-1)+t(i))- di(i)*(t(i)-taug(2*(i-1))))/(2*(taug(2*(i-1))-t(i-1))*(t(i)-t(i-1)));
		coeff(2*i-3,2) = di(i-1)-2*t(i-1)*coeff(2*i-3,1);
		coeff(2*i-3,3) = w1(i-1)-di(i-1)*t(i-1)+coeff(2*i-3,1)*t(i-1)^2;
		coeff(2*i-2,1) = (di(i)*(2*t(i)-t(i-1)-taug(2*(i-1)))-2*(w1(i)-w1(i-1))+di(i-1)*(taug(2*(i-1))-t(i-1)))/(2*(t(i)-taug(2*(i-1)))*(t(i)-t(i-1)));
		coeff(2*i-2,2) = di(i)-2*coeff(2*i-2,1)*t(i);
		coeff(2*i-2,3) = w1(i-1)+(taug(2*(i-1))-t(i-1))*(di(i-1)+(taug(2*(i-1))-t(i-1))*coeff(2*i-3,1))-taug(2*(i-1))*(coeff(2*i-2,2)+taug(2*(i-1))*coeff(2*i-2,1));
	end
end

clear functions
ts=max(abs(w1));  
%------------------------------------------------------------------
%SECOND PASS
%------------------------------------------------------------------
smurf = 10;
%function val2(k, coeff, taug) gives values for the value function and its
%derivatives.
fid=fopen('val2.m', 'w');
fprintf(fid, 'function v= val2(k, A, a)\n');
fprintf(fid, '%defines value and gradient of W given k and coefficient matrices\n');
fprintf(fid, 'if k<=a(1)\n');
fprintf(fid, 'j=1;\n');
fprintf(fid, 'else\n');
fprintf(fid, 'j=min(max(find(a<=k)), 44);\n');
fprintf(fid, 'end\n');
fprintf(fid, 'v1 = A(j,1)*k^2 + A(j, 2)*k + A(j,3);\n');
fprintf(fid, 'v2 = 2*A(j,1)*k + A(j, 2);\n');
fprintf(fid, 'v=[v1; v2];\n');
fclose(fid);

%vh2 gives high-state consumption that solves vh1(gamma), also the value of
%vh1(gamma) and its derivative.
fid=fopen('vh2.m', 'w');
fprintf(fid, 'function v=vh2(gam, A,a)\n');
fprintf(fid, 'global s p d h l\n');
fprintf(fid, 'if s>0&&gam<=d*a(45)\n');
fprintf(fid, 'Gone = gam/d;\n');
fprintf(fid, 'else\n');
fprintf(fid, 'Gone = a(45);\n');
fprintf(fid, 'end\n');
fprintf(fid, 'Gsol =0;\n');
fprintf(fid, 'if vhfoc(Gone, gam, A, a)<=0\n');
fprintf(fid, 'Gsol = Gone;\n');
fprintf(fid, 'else\n');
fprintf(fid, 'if gam <(1/s)*(.5+h)^s+d*a(1)\n');
fprintf(fid, 'Gtwo = a(1);\n');        
fprintf(fid, 'else\n');
fprintf(fid, 'Gtwo = (gam - (1/s)*(.5+h)^s)/d;\n');
fprintf(fid, 'end\n');
fprintf(fid, 'if vhfoc(Gtwo, gam, A, a)>=0\n');
fprintf(fid, 'Gsol = Gtwo;\n');
fprintf(fid, 'end\n');
fprintf(fid, 'end\n');
fprintf(fid, 'if Gsol ==0\n');
fprintf(fid, 'options = optimset(''LargeScale'', ''off'');\n');
fprintf(fid, 'Gsol = fzero(@vhfoc, [Gone Gtwo], options, gam, A, a);\n');
fprintf(fid, 'end\n');
fprintf(fid, 'x = (s*gam-s*d*Gsol)^(1/s);\n');
fprintf(fid, 'if x >= h+.5\n');
fprintf(fid, 'error(''dude 1 eating too much'')\n');
fprintf(fid, 'end\n');
fprintf(fid, 'fval = (1/s)*(1+h-x)^s+ [d 0]*val2(Gsol, A, a);\n');
fprintf(fid, 'der = ((1+h-x)/x)^(s-1);\n');
fprintf(fid, 'v = [x; fval; der];\n');
fclose(fid);
%auxiliary function vhfoc is the foc for use in vh2
fid=fopen('vhfoc.m', 'w');
fprintf(fid, 'function q = vhfoc(G, gam, A, a)\n');
fprintf(fid, 'global s p d h l\n');
fprintf(fid, '%for use in vh2\n');
fprintf(fid, 'q =-[0 1]*val2(G, A,a)-((s*gam-s*d*G)^(1/s)/(1+h-(s*gam-s*d*G)^(1/s)))^(1-s);\n');
fclose(fid);

% vlfoc is the first order condition to the problem of maximizing 1's utility in state L
% while holding 2 to autarky utility (needed for finding bhigh, Bhigh).
fopen('vlfoc.m', 'w');
fprintf(fid, 'function f=vlfoc(B, coeff, taug)\n');
fprintf(fid, 'global s p d h l\n');
fprintf(fid, 'b=1+l-(1/(1-d)-[d*s 0]*val2(B,coeff,taug))^(1/s);\n');
fprintf(fid, 'f=-[0 1]*val2(B, coeff, taug)-((1+l-b)/b)^(s-1);\n');
fclose(fid);
%as we know Bhigh<khigh, we only need to find the appropriate lower bound
%for that maximization problem. This is either klow or the B>klow that gives
%agent 2 his autarky utility when he has zero consumption in state L. The
%function Bbound gives this B when set to zero.
fopen('Bbound.m', 'w');
fprintf(fid, 'function f=Bbound(B, coeff, taug)\n');
fprintf(fid, 'global s p d h l\n');
fprintf(fid, 'f=-1/(s*(1-d))+[d 0]*val2(B, coeff, taug)+eps^s/s;\n');
fclose(fid);

%foctest2 is the first order condition to W(k).
fid=fopen('foctest2.m', 'w');
fprintf(fid, 'function f = foctest2(b, k,coeff,taug);\n');
fprintf(fid, 'global s p d h l\n');
fprintf(fid, 'f = p*[0 0 1]*vh2((k+(1-p)*((b+h-l)^s-b^s)/s),coeff,taug)*((b^(s-1)-(b+h-l)^(s-1))/(p*(b+h-l)^(s-1)+(1-p)*b^(s-1)))-(1+l-b)^(s-1)*(1/(p*(b+h-l)^(s-1)+(1-p)*b^(s-1)))-[0 1]*val2((k-(p*(b+h-l)^s+(1-p)*b^s)/s)/d,coeff,taug);\n');
fclose(fid);

% as in the first pass, for any b and k
%at which the revelation constraint holds, p2test1 shows whether P2L is
%slack, tight or violated.
fid=fopen('p2test2.m', 'w');
fprintf(fid, 'function f = p2test2(b, k,coeff,taug);\n');
fprintf(fid, 'global s p d h l xbar\n');
fprintf(fid, 'f = (1/s)*(1+l-b)^s+d*[1 0]*val2((k-(p*(b+h-l)^s+(1-p)*b^s)/s)/d,coeff,taug)-1/(s*(1-d));\n');
fclose(fid);

%As in the first pass, Blow and Bup (when set to zero) give the values of b
%for any k when B is set to klow/khigh.
fid=fopen('Blow.m', 'w');
fprintf(fid, 'function f = Blow(b, k,klow);\n');
fprintf(fid, 'global s p d h l\n');
fprintf(fid, 'f =k-d*klow-(p*(b+h-l)^s+(1-p)*b^s)/s;\n');
fclose(fid);
fid=fopen('Bup.m', 'w');
fprintf(fid, 'function f = Bup(b, k, khigh0);\n');
fprintf(fid, 'global s p d h l\n');
fprintf(fid, 'f =k-d*khigh0-(p*(b+h-l)^s+(1-p)*b^s)/s;\n');
fclose(fid);

while smurf >ts*(10^-12)
n=n+1
%reset values found in previous iteration to be starting values of next
%iteration

khigh0=khigh1;
w0=w1;
% in order to find the new upper bound khigh1, we need
% to find Bhigh, bhigh by maximizing 1's utility in state L
% while holding 2 to autarky utility. Therefore we need another foc,
% vlfoc.
Bone=khigh0;
if Bbound(klow, coeff, taug)>0
    Btwo = fzero ('Bbound', [klow, khigh0], options, coeff, taug);
else
    Btwo = klow;
end

if vlfoc(Btwo, coeff, taug)>=0
    Bhigh = Btwo;
elseif vlfoc(Bone, coeff, taug)<=0
    Bhigh = Bone;
else
    Bhigh = fzero('vlfoc', [Bone Btwo], options, coeff, taug);
end
bhigh = 1+l-((1/(1-d))-d*s*[1 0]*val2(Bhigh, coeff, taug))^(1/s);

% to test that the value function does not dip below the goal, we need
% to know that the revelation constraint holds at the new khigh

if (h^s)/s - ((bhigh + h-l)^s)/s +d*(khigh0 - Bhigh) < 0
	error ('RH violated at khigh');
end

kr = (1-p)*(bhigh^s)/s + p*((bhigh+h-l)^s)/s + d*Bhigh;
khigh1 = p*((h^s)/s  + d*khigh0) +(1-p)*((bhigh^s)/s + d*Bhigh);
betahigh = (bhigh^s)/s + d*Bhigh;

% we now evaluate the new value function at 20 points; note that to find
% the proper curvature, points are more closely spaced at the rim
% also note that for the chosen parameters, Bhigh>klow (we therefore omit
% programming for the case Bhigh=klow).
t(1) = klow;
hsol=vh2(((1-d*(1-p))*(h^s)+d*(1-p)*l^s)/(s*(1-d)),coeff,taug);
w1(1) = p*[0 1 0]*hsol+(1-p)*(1/s + d*[1 0]*val2(klow, coeff,taug));
bsol(1)=l;
q=1/(1+((1-p)/p)*(l/h)^(s-1));
di(1) = -q*[0 0 1]*hsol-(1-q)*(1/l)^(s-1);

t(22) = khigh1;
w1(23) = 1/(s*(1-d));
di(23) = -1/h^(s-1);
bsol(23) = bhigh;
    
for i = 2:11
	t(i) = klow + ((khigh1 - klow)/200)*(i - 1)^2;
end
for i=12:21
    t(i) = khigh1 - ((khigh1 - klow)/242)*(22-i)^2;
end
%specifying the solution for k>=kr
nr = min(find(t>= kr));
t(nr+1:23)=t(nr:22);
t(nr)=kr;

for i = nr:22
	hsol = vh2((t(i) -(1-p)*betahigh)/p,coeff,taug);
	w1(i) = p*[0 1 0]*hsol + (1-p)/(s*(1-d));
	di(i) = - [0 0 1]*hsol;
    bsol(i) = bhigh;
end

% next we check whether P2 binds at some k. For this we need ot identify
% for the chosen k the b that makes P2 exactly bind (if it exists), and 
% check whether the foc is greater or equal to zero. If k is
% so small that P2 is slack even at b=1, P2 will be slack at all smaller k.
% Otherwise, we know that P2 is slack at k<kr and b=bhigh, so that 
% the b that makes P2 exactly bind must exist.
% We also identify np=max(t<=kp).

for i = nr-1:(-1):2
    if p2test2(1+l-eps,t(i),coeff,taug)>=0
        np=i;
        break
    else
        k=t(i);
        bsol(i) = fzero('p2test2', [bhigh 1+l-eps], options, k,coeff,taug);
    end
    if foctest2(bsol(i),t(i),coeff,taug) > 0
		hsol = vh2(t(i)+(1-p)*((bsol(i)+h-l)^s-bsol(i)^s)/s,coeff,taug);
        w1(i) = p*[0 1 0]*hsol+(1-p)/(s*(1-d));
        vB=-[0 1]*val2((t(i)-(p*(bsol(i)+h-l)^s+(1-p)*bsol(i)^s)/s)/d,coeff,taug);
		di(i) = -p*[0 0 1]*hsol*(1+((1-p)*(bsol(i)^(s-1)-(bsol(i)+h-l)^(s-1)))/((1+l-bsol(i))^(s-1)/vB-p*(bsol(i)+h-l)^(s-1)-(1-p)*bsol(i)^(s-1)));
	else
		np = i;
		break
	end
end

%We proceed using foctest to find the solution for k<kp. We first identify
%the appropriate upper and lower bounds for b, bone and btwo.
%Having this solution we check that indeed B>=klow; if not, we set B=klow.
%We also identify nb as t(nb) is the largest element of t with B(k)=klow.

nb=1;
for i=2:np
    k=t(i);
    mick=0;
    if t(i)>= (p*(1+h-eps)^s+(1-p)*(1+l-eps)^s)/s+d*klow
        bone=1+l-eps;
    else
        bone=fzero('Blow', [l 1+l-eps], options, k,klow);
        mick=1;
    end
    if t(i)<=(p*h^s+(1-p)*l^s)/s+d*khigh0
        btwo=l;
    else
        btwo=fzero('Bup', [l 1+l-eps], options, k, khigh0);
    end
    if foctest2(bone,k,coeff,taug)>=0
        bsol(i)=bone;
    elseif foctest2(btwo,k,coeff,taug)<=0
        bsol(i)=btwo;
        mick=0;
    else
        bsol(i)=fzero('foctest2', [bone, btwo], options, k,coeff,taug);
        mick=0;
    end
    hsol=vh2(t(i)+(1-p)*((bsol(i)+h-l)^s-bsol(i)^s)/s,coeff,taug);
    if mick==1
        w1(i) = p*[0 1 0]*hsol+(1-p)*((1/s)*(1+l-bsol(i))^s+d*[1 0]*val2(klow,coeff,taug));
        nb=i;
    else
        w1(i) = p*[0 1 0]*hsol+(1-p)*((1/s)*(1+l-bsol(i))^s+d*[1 0]*val2((t(i)-(p*(bsol(i)+h-l)^s+(1-p)*bsol(i)^s)/s)/d,coeff,taug));
    end
    q=1/(1+((1-p)/p)*(bsol(i)/(bsol(i)+h-l))^(s-1));
    di(i) = -q*[0 0 1]*hsol-(1-q)*((1+l-bsol(i))/bsol(i))^(s-1);
end

clear functions

% we now compute the spline
for i=1:23
    taug(2*i-1)=t(i);
end
for i = 2:23
	if di(i - 1) < di(i)
		error ('value function not concave');
	end
end
for i = 2:23
    z(i) = (w1(i)-w1(i-1))/(t(i)-t(i-1));
    if z(i-1)<z(i)
        error('value function not concave at all');
    end
end
	
for i = 2:23	
	if .5*(di(i) + di(i-1)) == z(i)
		taug(2*(i-1)) = (t(i) + t(i-1))/2;
		coeff(2*i-3,1)=(di(i)-di(i-1))/(2*(t(i)-t(i-1)));
        coeff(2*i-3,2)= di(i-1)-2*t(i-1)*coeff(2*i-3,1);
        coeff(2*i-3,3) =w1(i-1)-di(i-1)*t(i-1)+t(i-1)^2;
		coeff(2*i-2,1)=coeff(2*i-3,1);
        coeff(2*i-2,2)=coeff(2*i-3,2);
        coeff(2*i-2,3)=coeff(2*i-3,3);
	else
		if (di(i-1)-z(i))*(di(i)-z(i))>=0
            error ('going weird again');
			taug(2*(i-1))=(t(i)+t(i-1))/2;
		elseif abs(di(i-1)-z(i))> abs(di(i)-z(i))
			taug(2*(i-1))=t(i-1)+(t(i)-t(i-1))*(di(i)-z(i))/(di(i)-di(i-1));
		else
			taug(2*(i-1))=t(i)+(t(i)-t(i-1))*(di(i-1)-z(i))/(di(i)-di(i-1));
		end
        coeff(2*i-3,1) = (2*(w1(i) - w1(i-1))-di(i-1)*(taug(2*(i-1))-2*t(i-1)+t(i))- di(i)*(t(i)-taug(2*(i-1))))/(2*(taug(2*(i-1))-t(i-1))*(t(i)-t(i-1)));
		coeff(2*i-3,2) = di(i-1)-2*t(i-1)*coeff(2*i-3,1);
		coeff(2*i-3,3) = w1(i-1)-di(i-1)*t(i-1)+coeff(2*i-3,1)*t(i-1)^2;
		coeff(2*i-2,1) = (di(i)*(2*t(i)-t(i-1)-taug(2*(i-1)))-2*(w1(i)-w1(i-1))+di(i-1)*(taug(2*(i-1))-t(i-1)))/(2*(t(i)-taug(2*(i-1)))*(t(i)-t(i-1)));
		coeff(2*i-2,2) = di(i)-2*coeff(2*i-2,1)*t(i);
		coeff(2*i-2,3) = w1(i-1)+(taug(2*(i-1))-t(i-1))*(di(i-1)+(taug(2*(i-1))-t(i-1))*coeff(2*i-3,1))-taug(2*(i-1))*(coeff(2*i-2,2)+taug(2*(i-1))*coeff(2*i-2,1));
	end
end
smurf = max(abs(w0-w1))
clear functions
end
%we now create policy functions for use in finding the stationary
%distribution.
%to find kB:
dud=1;
kone=t(nb);
ktwo=t(nb+1);
while dud>10^(-8)
    k=.5*(kone+ktwo);
    bone=fzero('Blow', [l 1+l-eps], options, k,klow);
    flag=foctest2(bone,k,coeff,taug);
    if flag<0
        ktwo=k;
    else
        kone=k;
    end
    dud=abs(flag)
end
kb=k;
bb=bone;
nb=nb+1;
%having kb, we can now create data vectors. 
for i=1:nb-1
    hsol=vh2(t(i)+(1-p)*((bsol(i)+h-l)^s-bsol(i)^s)/s,coeff,taug);
    dat(i,1)=(t(i)+(1-p)*((bsol(i)+h-l)^s-bsol(i)^s)/s-(hsol(1))^s/s)/d;
    dat(i,2)=klow;
    int(i)=t(i);
end
hsol=vh2(kb+(1-p)*((bb+h-l)^s-bb^s)/s,coeff,taug);
dat(nb,1)=(kb+(1-p)*((bb+h-l)^s-bb^s)/s-(hsol(1))^s/s)/d;
dat(nb,2)=klow;
int(nb)=kb;

for i=nb:nr-1
    hsol=vh2(t(i)+(1-p)*((bsol(i)+h-l)^s-bsol(i)^s)/s,coeff,taug);
    dat(i+1,1)=(t(i)+(1-p)*((bsol(i)+h-l)^s-bsol(i)^s)/s-(hsol(1))^s/s)/d;
    dat(i+1,2)=(t(i)-(p*(bsol(i)+h-l)^s+(1-p)*bsol(i)^s)/s)/d;
    int(i+1)=t(i);
end
dat(nr+1,1)=kr;
dat(nr+1,2)=Bhigh;
int(nr+1)=kr;

%normalize values of kappa
dat=(dat-klow)*(1/(khigh0-klow));
int=(int-klow)*(1/(khigh0-klow));
Bhigh=(Bhigh-klow)*(1/(khigh0-klow));
kb=(kb-klow)*(1/(khigh0-klow));
kr=(kr-klow)*(1/(khigh0-klow));

%give a visual check on derived policy functions
plot(int, dat)
%we now produce the approximate ergodic set and the corresponding
%stationary distribution.  
grid=[0:1/20000:kr kr]';
lra=length(grid);
pig=max(find(grid<=kb));
test(1:pig)=interp1(int(1:nb),dat(1:nb,1),grid(1:pig),'linear');
test(pig+1:lra)=interp1(int(nb:nr+1),dat(nb:nr+1,1),grid(pig+1:lra),'linear');
if max(test(1:lra-2)'-grid(1:lra-2))<1/20000
    error('grid too coarse')
end
trans=zeros(lra-1,1);
test=0;
erg1=0;
n=0;
%we first need to populate K, the set of attainable utilities, so that we
%have an approximation to the ergodic set that spans K.
while test<1 
    n=n+1
    erg0=erg1;
    lmu=length(erg0);
    if lmu>=2^20
        C=erg0(lmu);
        [c,pc]=sort([grid;erg0(1:lmu-1)]);
        gl=find(pc<=lra);
        for i=1:lra-1
            if gl(i+1)==gl(i)+1
                trans(i)=2;
            else
                trans(i)=erg0(pc(gl(i)+1)-lra,1);
            end
        end
        gl=find(trans<2);
        erg0=trans(gl);
        erg0(end+1)=C;
        lmu=length(erg0);
    end
    erg1=repmat(erg0,2,1);
    pig=max(find(erg0<=kb));
    erg1(1:pig)=zeros(pig,1);
    erg1(pig+1:2*pig)=interp1(int(1:nb),dat(1:nb,1),erg0(1:pig),'linear');
    erg1(2*pig+1:pig+lmu)=interp1(int(nb:nr+1),dat(nb:nr+1,2),erg0(pig+1:lmu),'linear');
    erg1(pig+lmu+1:2*lmu)=interp1(int(nb:nr+1),dat(nb:nr+1,1),erg0(pig+1:lmu),'linear');
    erg1=sort(erg1);
	gl=find(diff(erg1)==0);
    erg1(gl+1)=[];
    max(erg1)
    if max(erg1)>=grid(end-1)&&length(erg1)>=100000
        test=1;
	end
end
%we then fix the size of the ergodic set (to be able to generate
%eigenvectors after - if K is too large, eigenvalue routines fail)
clear erg0 
C=erg1(end);
[c,pc]=sort([grid;erg1(1:end-1)]);
gl=find(pc<=lra);
for i=1:lra-1
    if gl(i+1)==gl(i)+1
        trans(i)=2;
    else
        trans(i)=erg1(pc(gl(i)+1)-lra,1);
    end
end
gl=find(trans<2);
erg1=trans(gl);
erg1(end+1)=C;
if erg1(end)<kr
    erg1(end+1)=kr;
end
%we now generate the transition matrix
clear trans grid;
lmu=length(erg1);
B=zeros(lmu,1);
G=zeros(lmu,1);
bi=zeros(1,lmu);
gi=zeros(1,lmu);
row=zeros(1,4*lmu);
col=zeros(1,4*lmu);
pr=zeros(1,4*lmu);
pig=max(find(erg1<=kb));
B(1:pig)=zeros(pig,1);
B(pig+1:lmu)=interp1(int(nb:nr+1),dat(nb:nr+1,2),erg1(pig+1:lmu),'linear');
G(1:pig)=interp1(int(1:nb),dat(1:nb,1),erg1(1:pig),'linear');
G(pig+1:lmu)=interp1(int(nb:nr+1),dat(nb:nr+1,1),erg1(pig+1:lmu),'linear');
[c, pc]=sort([erg1;B]);
gl=find(pc<=lmu);
for i=1:lmu-1
    if gl(i+1)>gl(i)+1
        x=pc(gl(i)+1:gl(i+1)-1)-lmu;
        bi(x)=i;
    end
end
[c, pc]=sort([G;erg1]);
gl=find(pc>lmu);
for i=2:lmu
    if gl(i)>gl(i-1)+1
        x=pc(gl(i-1)+1:gl(i)-1);
        gi(x)=i;
    end
end
for i=1:lmu
    i
    col(4*i-3:4*i)=i;
    row(4*i-3:4*i)=[bi(i) bi(i)+1 gi(i)-1 gi(i)];
    pr(4*i-3:4*i-2)=[erg1(bi(i)+1)-B(i) B(i)-erg1(bi(i))]*(1-p)/(erg1(bi(i)+1)-erg1(bi(i)));
	pr(4*i-1:4*i)=[erg1(gi(i))-G(i) G(i)-erg1(gi(i)-1)]*p/(erg1(gi(i))-erg1(gi(i)-1));
	pr(4*i-3:4*i)=pr(4*i-3:4*i)/sum(pr(4*i-3:4*i));
    if bi(i)+1==gi(i)-1 && pr(4*i-2)>0 && pr(4*i-1)>0
        pr(4*i-3:4*i)=[.5 0 0 .5];
    end
end
gl=find(pr==0);
pr(gl)=[];
row(gl)=[];
col(gl)=[];
clear pc c gl B G bi gi;
TP=sparse(row,col,pr,lmu,lmu);
if min(sum(TP))<1
    error('columns of TP don''t add up to 1')
end
clear row col pr;
[mu,di]=eigs(TP,1,1-1e-6);
if mu(1)<0
    mu=-mu;
end
if min(mu)<0
    error('eigenvector negative');
end
mu=mu*(1/sum(mu));
plot(erg1,mu)