%SOLVEX, a script file (not a function).
%
%   SOLVEX solves the Gaussian exponential quadratic
%   social planning problem, 
%
%     V(x(t)) = max { -.5  [(s(t)-b(t)).(s(t)-b(t)) + l(t)^2] 
%        
%           + beta *(2/sig) * log E(t) exp  ( sig * V(x(t+1)/2 )}
%
%   subject to
%
%      phic * c(t) + phig * g(t) + phii * i(t) = gamma * k(t-1) + d(t)
%                              g(t).g(t) = l(t)^2
%                                   k(t) = deltak * k(t-1) + thetak * i(t)
%                                   h(t) = deltah * h(t-1) + thetah * c(t)
%                                   s(t) = lambda * h(t-1) + pih * c(t)
%                                 z(t+1) = a22 * z(t) + c2 * w(t+1)
%                                   b(t) = ub * z(t), d(t) = ud * z(t)
%   Where x(t) = [h(t-1)',k(t-1)',z(t)']', the solution of the problem is
%
%        x(t+1) = ao * x(t) + c * w(t+1)
%        j(t) = sj * x(t)
%   where j = k, h, c, s, g, i, b, and d; w(t+1) is a Gaussian
%   white noise with unit contemporaneous covariance matrix.
%   The program also computes
%   Lagrange multipliers
%
%        muj = mj * x(t)
%   for j=k,c,h,s, and i.
%   When sig = 0, this collapses to linear quadratic; when sig is LESS
%   than zero, we have increased risk-aversion relative to linear
%   quadratic.

%STEP 1:
%Computation of the dimension of the structural parameter matrices that the
%user has specified.  The user is told the dimension of the various
%matrices in order to help the user verify that the matrices are inputted
%correctly.

[nb,nh]=size(lambda);
[nd,nc]=size(phic);
[nz,nw]=size(c2);
[junk,ng]=size(phig);
[nk,ni]=size(thetak);


%STEP 2:
%Creation of the composite state vector x(t)'=[h(t-1)',k(t-1)',z(t)'] and its
%law of motion:  x(t+1) = A*x(t) + B*i(t) + C*w(t+1).

  phiin = inv([phic phig]);
  phiinc = [eye(nc) zeros(nc,ng)]*phiin;
  phiing = [zeros(ng,nc) eye(ng)]*phiin;
  b11 = - thetah*phiinc*phii;
  a1 = thetah*phiinc*gamma;
  a12 = [thetah*phiinc*ud; zeros(nk,nz)];

    %Creation of the A matrix

     a11=[deltah a1 ; zeros(nk,nh) deltak];
     a=[a11 a12; [zeros(nz,nk+nh) a22]];

    %Creation of the B matrix

     b1=[b11; thetak];
     b=[b1; zeros(nz,ni)];

    %Creation of the C matrix

     c = [zeros(nk+nh,nw);c2];


%Scaling of the law of motion by the square root of discount factor to account
%for the discounting in the optimization problem

   a=a*sqrt(beta);

   a11=a11*sqrt(beta);

   b1=b1*sqrt(beta);

   b=b*sqrt(beta);



%STEP 3:
%Creation of the matrices Q and R for the objective function of the optimal
%linear regulator problem.  The matrices Q, R, S are first created where
%the matrix S allows for the presence of cross products between the control
%and the state in the objective function.  The control is then transformed in
%the last step to remove the cross product and to put the problem in the proper
%form of the optimal linear regulator problem.


%STEP 3(i):
%Form matrices Q, R, S where S contains cross-products

   n=ni+nh+nk+nz;
   nx = nh+nk+nz;

   g1 = phiing*[-phii zeros(nd,nh) gamma ud];
   g2 = [-pih*phiinc*phii lambda pih*phiinc*gamma -ub+pih*phiinc*ud];

   r=g1'*g1 + g2'*g2;

   q=r(1:ni,1:ni);

   s=r(1:ni,1+ni:n);

   r=r(1+ni:n,1+ni:n);


%STEP 3(ii)
%Transform control to remove cross-product

   %First check if q is nonsingular

if rank(q) < max(size(q));

   error('EXECUTION STOPPED.  The 1-1 block of the R matrix is singular')

   end

   qi=inv(q);

   a=a-b*qi*s;

   r=r-s'*qi*s;

   a11=a(1:nh+nk,1:nh+nk);

   r1=r(1:nh+nk,1:nh+nk);


%STEP 4:
%Solve the Riccati equation useing two doubling algorithms.                  
%We exploit the fact that the solution to the Riccati equation can be solved
%in parts.  


   disp('Calculating, please wait')

%STEP 4(i):
%Calculate the 11 block of the Riccati equation.  The result is a matrix kg
%such that kg'*[h(t), k(t)] is the feedback part of the control law.             


   [k,v,vt]=doublex(a',b',r,q,c,-sig);

%    remember that doublex.m solves a MIN problem, which requires
%    us to reverse the sign convention for sig

%   a11o=a11-b1*kg';      %a11o gives the 11 (all endogenous variables)
                         %part of ao.                                          

%STEP 4(ii):
%Calculate the 12 block of the Riccati equation.  The results is the matrix
%kg2 such that kg2*z(t) is the feedforward part of the control law.              

%   a22t=a(nk+nh+1:nx,nk+nh+1:nx);

%   r2=r(1:nh+nk,nh+nk+1:nx);

%   a12=a(1:nh+nk,nh+nk+1:nx);

%   qopt=inv(q+b1'*v11*b1);

 %  r12=r2-a11'*v11*b1*qopt*b1'*v11*a12+a11'*v11*a12;

%   v12=doublej2(a11o',r12,a22t',eye(nz));

%   kg2=qopt*b1'*(v12*a22t+v11*a12);

%STEP 4(iii):
%Calculate the matrix ao which gives the reduced form law of motion for
%x(t) of the form:  x(t+1) = ao*x(t) + c*w(t+1).  The matrix ao also has
%the square root of beta rescalling unwound.


   ao=(a-b*k')/sqrt(beta); %Square root of beta
                                    %unwound.



   a11o=ao(1:nk+nh,1:nk+nh)          %Square root of beta unwound in the
                                    %endogenous state variable part of ao


%STEP 5:
%Display eigenvalue calculations for the user.  

%Endogenous eigenvalues.                                                         


   endo=eig(a11o);


%Exogenous eigenvalues.

   exo=eig(a22);


%STEP 6:
%Calculation of the solution matrices.                                      

    sh = ao(1:nh,1:nx);
    sk = ao(nh+1:nh+nk,1:nx);
    si = -k' - qi*s;
    skl = [zeros(nk,nh) eye(nk,nk) zeros(nk,nz)];
    sd = [zeros(nd,nh+nk),ud];
    sb = [zeros(nb,nh+nk),ub];
    sc = phiinc*(-phii*si + gamma*skl + sd);
    sg = phiing*(-phii*si + gamma*skl + sd);
    ss = lambda*[eye(nh,nh) zeros(nh,nk+nz)] + pih*sc;


disp('The matrix ao has been calculated for the law of ')
disp('motion of the entire state vector.  This matrix satisfies')
disp('        x(t+1) = ao*x(t) + c*w(t+1). ')
disp('The endogenous eigenvalues are in the vector endo,')
disp('and the exogenous eigenvalues are in the vector exog.')
disp('   ')
disp('   ')

disp('The solution to the model is given by c(t) = sc*x(t),')
disp(' g(t) = sg*x(t), h(t) = sh*x(t), i(t) = si*x(t)')
disp(' k(t) = sk*x(t), and s(t) = ss*x(t).')
disp('The matrices sc, sg, sh, si, sk, and ss have now been computed')
disp('and can be used in other matlab programs.')
disp('The matrices sb and sd are constructed so that b(t) = sb*x(t)')
disp('and d(t) = sd*x(t) and can be used in other matlab programs.')
%Calculating the Lagrange Multipliers%


v1=vt(1:nh+nk,1:nx);

    mh = - beta*[eye(nh) zeros(nh,nk)]*[v1]*ao;
    mk = - beta*[zeros(nk,nh) eye(nk)]*[v1]*ao;
    ms = sb - ss;
    mc = thetah'*mh + pih'*ms;
    mg = sg;
    mi = thetak'*mk;
    md = phiin'*[mc; -sg];
										
										
disp('     ')
disp('     ')

disp('The shadow price vectors satisfy Mc(t) = mc*x(t), Mg(t) = mg*x(t),')
disp('Mh(t) = mh*x(t), Mi(t) = mi*x(t), Mk(t) = mk*x(t), Ms(t) = ms*x(t),')
disp('and Md(t) = md*x(t).  The matrices of these linear combinations can')
disp('be used in other matlab programs.')
disp('     ')
disp('     ')

disp('Your equilibrium has been calculated')
disp('     ')
