% SOLVES.M
%   Solves.m is a script file that solves a periodic version of the
%   social planning problem, to maximize
%
%   -.5 SUM [(s(t)-b(t)).(s(t)-b(t)) + l(t)^2]
%   subject to
%   phics(t)*c(t)+phigs(t)*g(t)+phiis(t)*i(t) = gammas(t)*k(t-1) +d(t)
%                   g(t).g(t) = l(t)^2
%                 k(t) = deltaks(t)*k(t-1) + thetaks(t)*i(t)
%                 h(t) = deltahs(t)*h(t-1) + thetahs(t)*c(t)
%                 s(t) = lambdas(t)*h(t-1) + pihs(t)*c(t)
%               z(t+1) = a22s(t)*z(t) + c2s(t) * w(t+1)
%                 b(t) = ub * z(t), d(t) = ud * z(t)
%   where s(t+p) = s(t) is a periodic function into the integers, and p is
%   the period of the model.
%   Where x(t) = [h(t)', k(t)', z(t)'], the solution of the problem is
%
%        x(t+1) = aos(t) * x(t) + cs(t) * w(t+1)
%          j(t) = sjs(t) * x(t)
%   where j = k,h,c,s,g, i, b, and d.  The program also computes
%   the matrices mjs(t) that pin down the Lagrange multipliers mujs(t) in
%         mujs(t) = mjs(t) * x(t)
%    for j = k,c,h,s, and i.
										
load seas4        % loads in shells of matrices
p=4;            % sets value of periodicity
seas1               % loads in a seasonal matrices

gamma1(1)=.13;
gamma4(1)=.08;										
disp('   ')
disp(' Calculating, please wait')
										
[nb,nh]=size(lambda1);
[nd,nc]=size(phic1);
[nz,nw]=size(c21);
[junk,ng]=size(phig1);
[nk,ni]=size(thetak1);
for j=1:p
n=int2str(j);
eval(['phiin',n,'=inv([phic',n,' phig',n,']);'])
eval(['phiinc',n,'=[eye(nc) zeros(nc,ng)]*phiin',n,';'])
eval(['phiing',n,'=[zeros(ng,nc) eye(ng)]*phiin',n,';'])
eval(['b11',n,'=-thetah',n,'*phiinc',n,'*phii',n,';'])
eval(['ai',n,'=thetah',n,'*phiinc',n,'*gamma',n,';'])
eval(['a12',n,'=[thetah',n,'*phiinc',n,'*ud;zeros(nk,nz)];'])
										
%  Create A matrices
										
eval(['a11',n,'=[deltah',n,' ai',n,';zeros(nk,nh) deltak',n,'];'])
eval(['A',n,'=[a11',n,' a12',n,';[zeros(nz,nk+nh) a22',n,']];'])
										
%  Create B matrices
eval(['BB1',n,'=[b11',n,';thetak',n,'];'])
eval(['B',n,'=[BB1',n,';zeros(nz,ni)];'])
										
%   Create C matrices
										
eval(['C',n,'=[zeros(nk+nh,nw);c2',n,'];'])
										
%   Scaling
										
eval(['A',n,'=sqrt(beta)*A',n,';'])
eval(['a11',n,'=sqrt(beta)*a11',n,';'])
eval(['BB1',n,'=sqrt(beta)*BB1',n,';'])
eval(['B',n,'=sqrt(beta)*B',n,';'])
										
										
nn=ni+nh+nk+nz;
nx=nh+nk+nz;
										
eval(['g1',n,'=phiing',n,'*[-phii',n,' zeros(nd,nh) gamma',n,' ud];'])
eval(['g2',n,'=[-pih',n,'*phiinc',n,'*phii',n,' lambda',n,' pih',n,...
      '*phiinc',n,'*gamma',n,' -ub+pih',n,'*phiinc',n,'*ud];'])
eval(['r',n,'=g1',n,'''*g1',n,'+g2',n,'''*g2',n,';'])
eval(['q',n,'=r',n,'(1:ni,1:ni);'])
eval(['s',n,'=r',n,'(1:ni,1+ni:nn);'])
eval(['r',n,'=r',n,'(1+ni:nn,1+ni:nn);'])
										
										
%  Transform to remove cross-products
										
eval(['qi',n,'=inv(q',n,');'])
eval(['A',n,'=A',n,'-B',n,'*qi',n,'*s',n,';'])
eval(['r',n,'=r',n,'-s',n,'''*qi',n,'*s',n,';'])
eval(['a11',n,'=A',n,'(1:nh+nk,1:nh+nk);'])
eval(['ri',n,'=r',n,'(1:nh+nk,1:nh+nk);'])
end
for j=1:p
n=int2str(j);
eval(['aa',n,'=a11',n,';'])
eval(['gg',n,'=ri',n,';'])
%eval(['xx=b1',n,';'])
%eval(['yy=q',n,';']);
%zz=xx';
%cc=xx*(yy\zz);
%eval(['b',n,'=cc;'])
%eval(['h',n,'=b1',n,'*(q',n,'\zz);'])
eval(['bb',n,'=BB1',n,'*(q',n,'\BB1',n,''');'])
end
										
%  Use symplectic matrices to solve dynamic programming problem
										
%  First create Z(p-1)*Z(p)
										
np=int2str(p);
pm1=int2str(p-1);
eval(['[AA,Bd,GG]=mult(aa',np,',bb',np,',gg',np,',aa',pm1,',bb',pm1,',gg'...
	 ,pm1,');'])
%  Now create Z(1)*Z(2)*...*Z(1)
										
for j=p-2:-1:1;
n=int2str(j);
eval(['[AA,Bd,GG]=mult(AA,Bd,GG,aa',n,',bb',n,',gg',n,');'])
end
										
%  Compute P1
										
[AA1,Bd1,P1]=double3(AA,Bd,GG);
										
%  Compute Pp, P(p-1),..., P2
										
eval(['[AA',np,',Bd',np,',P',np,']=mult(AA1,Bd1,P1,aa',np,',bb',np,',gg',np,...
    ');'])
for j=p-1:-1:2;
n1=int2str(j+1);
n=int2str(j);
eval(['[AA',n,',Bd',n,',P',n,']=mult(AA',n1,',Bd',n1,',P',n1,',aa',n,',bb',n,..
    ',gg',n,');'])
end
np1=int2str(p+1);
eval(['P',np1,'=P1;'])       %  Set P(p+1) = P1
for j=1:p;
n=int2str(j);
np1=int2str(j+1);
eval(['qopt',n,'=inv(q',n,'+BB1',n,'''*P',np1,'*BB1',n,');'])
eval(['f',n,'=qopt',n,'*BB1',n,'''*P',np1,'*a11',n,';'])
end
										
										
										
										
										
%  Calculate closed loop matrices
for j=1:p
n=int2str(j);
eval(['a11o',n,'=a11',n,'-BB1',n,'*f',n,';'])
end
										
% calculate (1,2) block of Riccati equation
for j=1:p
n=int2str(j);
eval(['a22t',n,'=A',n,'(nk+nh+1:nx,nk+nh+1:nx);'])
eval(['r2',n,'=r',n,'(1:nh+nk,nh+nk+1:nx);'])
eval(['a12',n,'=A',n,'(1:nh+nk,nh+nk+1:nx);'])
end
										
										
										
npp=int2str(p+1);
eval(['P',npp,'=P1;'])  % sets value function periodically
for j=1:p
n=int2str(j);
np1=int2str(j+1);
eval(['D',n,'=r2',n,'+a11',n,'''*P',np1,'*a12',n,...
     '-a11',n,'*P',np1,'*BB1',n,'*qopt',n,'*BB1',n,'''*P',np1,...
    '*a12',n,';'])
eval(['G',n,'=a11',n,'-BB1',n,'*qopt',n,'*BB1',n,'''*P',np1,'*a11',n,';'])
eval(['H',n,'=a22t',n,';'])
end
										
%  Create H,G,D for seasonal doubling of v12
p=p;
h1=H1;
g1=G1;
npp=int2str(p);
eval(['d1=D',npp,';'])
for j=2:p
n=int2str(j);
npj=int2str(p+1-j);
nm1=int2str(j-1);
eval(['g',n,'=G',n,'*g',nm1,';'])
eval(['h',n,'=H',n,'*h',nm1,';'])
eval(['d',n,'=G',npj,'''*d',nm1,'*H',npj,'+D',npj,';'])
end
eval(['G=g',n,';'])
eval(['H=h',n,';'])
eval(['D=d',n,';'])
										
v121=doublej2(G',D,H',eye(nz));
										
pp1=int2str(p+1);
eval(['v12',pp1,'=v121;'])
for j=p:-1:2
n=int2str(j);
np1=int2str(j+1);
eval(['v12',n,'=D',n,'+G',n,'''*v12',np1,'*H',n,';'])
end
										
										
for j=1:p
n=int2str(j);
np1=int2str(j+1);
eval(['f2',n,'=qopt',n,'*BB1',n,'''*(v12',np1,'*a22t',n,'+P',np1,...
      '*a12',n,');'])
end
										
										
%   Create closed loop matrices
%
for j=1:p
n=int2str(j);
eval(['ao',n,'=(A',n,'-B',n,'*[f',n,' f2',n,'])/sqrt(beta);'])
eval(['a11o',n,'=a11o',n,'/sqrt(beta);'])
end
										
										
%  Create S's giving variables as linear functions of state
										
for j=1:p
n=int2str(j);
eval(['sh',n,'=ao',n,'(1:nh,1:nx);'])
eval(['sk',n,'=ao',n,'(nh+1:nh+nk,1:nx);'])
eval(['si',n,'=-[f',n,' f2',n,']-qi',n,'*s',n,';'])
eval(['sk1',n,'=[zeros(nk,nh) eye(nk,nk) zeros(nk,nz)];'])
eval(['sd',n,'=[zeros(nd,nh+nk),ud];'])
eval(['sb',n,'=[zeros(nb,nh+nk),ub];'])
eval(['sc',n,'=phiinc',n,'*(-phii',n,'*si',n,...
     '+gamma',n,'*sk1',n,'+sd',n,');'])
eval(['sg',n,'=phiing',n,'*(-phii',n,'*si',n,...
     '+gamma',n,'*sk1',n,'+sd',n,');'])
eval(['ss',n,'=lambda',n,'*[eye(nh,nh) zeros(nh,nk+nz)]'...
     '+pih',n,'*sc',n,';'])
end
										
										
%  create M's for Lagrange multipliers
for j=1:p
n=int2str(j);
np1=int2str(j+1);
eval(['mh',n,'=-beta*[eye(nh) zeros(nh,nk)]*[P',np1,...
     ' v12',np1,']*ao',n,';'])
eval(['mk',n,'=-beta*[zeros(nk,nh) eye(nk)]*[P',np1,...
     ' v12',np1,']*ao',n,';'])
eval(['ms',n,'=sb',n,'-ss',n,';'])
eval(['mc',n,'=thetah',n,'''*mh',n,'+pih',n,'*ms',n,';'])
eval(['mg',n,'=sg',n,';'])
eval(['mi',n,'=thetak',n,'''*mk',n,';'])
eval(['md',n,'=phiin',n,'*[mc',n,'; -sg',n,'];'])
end
disp('    ')
disp(' Your equilibrium has been calculated')
disp('   ')
disp('The shadow price vectors satisfy Mc(t) = mcs(t)*x(t), Mg(t) = ')
disp('  mgs(t) * x(t), Mh(t) = mhs(t)*x(t), Mi(t) = mis(t)*x(t), ')
disp('  Mk(t) = mks(t) * x(t), Ms(t) = mss(t) * x(t), and ')
disp('  Mh(t) = mhs(t) * x(t).  These matrices are stored and can be ')
disp('  used in other Matlab programs.')
disp('  Similarly, c(t) = Scs(t) * x(t), and so on.  The matrices')
disp('  Sjs(t) are stored.')
										
