% Main file to generate intial guess for solution to model A2
% (symmetric CRRA utility function over cons & labor, symmetric CD production function)
% Approximation is done in levels (k) and logs (a)

clear all

cpu0=cputime;

global N beta alpha delta sigma rho phi A gamma tau eta b
global kmin kmax amin amax
global xgrid
global k1_mat k2_mat k3_mat k4_mat k5_mat k6_mat k7_mat k8_mat K_mat
global a1_mat a2_mat a3_mat a4_mat a5_mat a6_mat a7_mat a8_mat
global ap1_mat ap2_mat ap3_mat ap4_mat ap5_mat ap6_mat ap7_mat ap8_mat
global y1_mat y2_mat y3_mat y4_mat y5_mat y6_mat y7_mat y8_mat Y_mat

%%%%%%%%%%%%%%%%%%%%%%
% 1. Getting started %
%%%%%%%%%%%%%%%%%%%%%%

approx=1;   % set order of approximation for perturbation solution (initilization)

% 1) Calibration
N=8;
beta=0.99;
alpha=0.36;
delta=0.025;
sigma=0.01;
rho=0.95;
phi=0.5;
A=(1-beta)/(alpha*beta);
% gammalow=0.25;
% gammahigh=1;
% gamma=[gammalow:(gammahigh-gammalow)/(N-1):gammahigh];
gamma=ones(1,N);

% etalow=0.1;
% etahigh=1;
% eta=-[etalow:(etahigh-etalow)/(N-1):etahigh];
eta=-ones(1,N);


% 2) Compute the steady state, the Negishi weights and the labor weight in
% the utility function

% a)
kss=1;
lss=1;
Kss=N*kss;
Lss=N*lss;
css=A*kss^alpha;
Css=N*css;
lnass=0;
ass=1;
thetass=1;
Yss=A*Kss^alpha*Lss^(1-alpha);

% b) Compute Negishi weights
tau=zeros(1,N);
for j=1:N
    tau(j)=css^(1/gamma(j));
end
muss=1;

% c) Labor weights
b=zeros(1,N);
for j=1:N
    b(j)=(1-alpha)*A^(1-1/gamma(j));
end


% 3) Compute the steady state in logs
lnkss=log(kss);
lnKss=log(Kss);
lncss=log(css);
lnCss=log(Css);
lnthetass=log(thetass);
lnYss=log(Yss);
lnLss=log(Lss);
lnlss=log(lss);
lnmuss=log(muss);



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 2. Construct first-order approximation using method by Klein (2000)  %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 1) Define symbolic variables
syms c1 lna1 k1 theta1 cp1 lnap1 kp1 thetap1 l1 lp1
syms c2 lna2 k2 theta2 cp2 lnap2 kp2 thetap2 l2 lp2
syms c3 lna3 k3 theta3 cp3 lnap3 kp3 thetap3 l3 lp3
syms c4 lna4 k4 theta4 cp4 lnap4 kp4 thetap4 l4 lp4
syms c5 lna5 k5 theta5 cp5 lnap5 kp5 thetap5 l5 lp5
syms c6 lna6 k6 theta6 cp6 lnap6 kp6 thetap6 l6 lp6 
syms c7 lna7 k7 theta7 cp7 lnap7 kp7 thetap7 l7 lp7 
syms c8 lna8 k8 theta8 cp8 lnap8 kp8 thetap8 l8 lp8 
syms C Cp K Kp Y Yp CAC CACp mu mup L Lp

% 1) Define symbolic variables
syms lnc1 lna1 lnk1 lntheta1 lncp1 lnap1 lnkp1 lnthetap1 lnl1 lnlp1
syms lnc2 lna2 lnk2 lntheta2 lncp2 lnap2 lnkp2 lnthetap2 lnl2 lnlp2
syms lnc3 lna3 lnk3 lntheta3 lncp3 lnap3 lnkp3 lnthetap3 lnl3 lnlp3
syms lnc4 lna4 lnk4 lntheta4 lncp4 lnap4 lnkp4 lnthetap4 lnl4 lnlp4
syms lnc5 lna5 lnk5 lntheta5 lncp5 lnap5 lnkp5 lnthetap5 lnl5 lnlp5 
syms lnc6 lna6 lnk6 lntheta6 lncp6 lnap6 lnkp6 lnthetap6 lnl6 lnlp6
syms lnc7 lna7 lnk7 lntheta7 lncp7 lnap7 lnkp7 lnthetap7 lnl7 lnlp7 
syms lnc8 lna8 lnk8 lntheta8 lncp8 lnap8 lnkp8 lnthetap8 lnl8 lnlp8 
syms lnC lnCp lnK lnKp lnY lnYp CAC CACp lnmu lnmup lnL lnLp

% 2) Construct the system of equilibrium conditions

f1=exp(lnc1)^(-1/gamma(1))*[1+phi*(exp(lnkp1)/exp(lnk1)-1)]-beta*exp(lncp1)^(-1/gamma(1))*[exp(lnthetap1)+exp(lnap1)*A*alpha*exp(lnkp1)^(alpha-1)*exp(lnlp1)^(1-alpha)];
f1b=b(1)*[exp(lnl1)^(alpha-1/eta(1))]-[exp(lnc1)^(-1/gamma(1))]*exp(lna1)*(1-alpha)*A*exp(lnk1)^alpha;
f2=exp(lntheta1)-[1+phi*[1+0.5*(exp(lnkp1)/exp(lnk1)-1)]*(exp(lnkp1)/exp(lnk1)-1)];
f3=tau(1)*exp(lnc1)^(-1/gamma(1))-exp(lnmu);
f4=lnap1-rho*lna1;

f5=exp(lnc2)^(-1/gamma(2))*[1+phi*(exp(lnkp2)/exp(lnk2)-1)]-beta*exp(lncp2)^(-1/gamma(2))*[exp(lnthetap2)+exp(lnap2)*A*alpha*exp(lnkp2)^(alpha-1)*exp(lnlp2)^(1-alpha)];
f5b=b(2)*[exp(lnl2)^(alpha-1/eta(2))]-[exp(lnc2)^(-1/gamma(2))]*exp(lna2)*(1-alpha)*A*exp(lnk2)^alpha;
f6=exp(lntheta2)-[1+phi*[1+0.5*(exp(lnkp2)/exp(lnk2)-1)]*(exp(lnkp2)/exp(lnk2)-1)];
f7=tau(2)*exp(lnc2)^(-1/gamma(2))-exp(lnmu);
f8=lnap2-rho*lna2;

f9=exp(lnc3)^(-1/gamma(3))*[1+phi*(exp(lnkp3)/exp(lnk3)-1)]-beta*exp(lncp3)^(-1/gamma(3))*[exp(lnthetap3)+exp(lnap3)*A*alpha*exp(lnkp3)^(alpha-1)*exp(lnlp3)^(1-alpha)];
f9b=b(3)*[exp(lnl3)^(alpha-1/eta(3))]-[exp(lnc3)^(-1/gamma(3))]*exp(lna3)*(1-alpha)*A*exp(lnk3)^alpha;
f10=exp(lntheta3)-[1+phi*[1+0.5*(exp(lnkp3)/exp(lnk3)-1)]*(exp(lnkp3)/exp(lnk3)-1)];
f11=tau(3)*exp(lnc3)^(-1/gamma(3))-exp(lnmu);
f12=lnap3-rho*lna3;

f13=exp(lnc4)^(-1/gamma(4))*[1+phi*(exp(lnkp4)/exp(lnk4)-1)]-beta*exp(lncp4)^(-1/gamma(4))*[exp(lnthetap4)+exp(lnap4)*A*alpha*exp(lnkp4)^(alpha-1)*exp(lnlp4)^(1-alpha)];
f13b=b(4)*[exp(lnl4)^(alpha-1/eta(4))]-[exp(lnc4)^(-1/gamma(4))]*exp(lna4)*(1-alpha)*A*exp(lnk4)^alpha;
f14=exp(lntheta4)-[1+phi*[1+0.5*(exp(lnkp4)/exp(lnk4)-1)]*(exp(lnkp4)/exp(lnk4)-1)];
f15=tau(4)*exp(lnc4)^(-1/gamma(4))-exp(lnmu);
f16=lnap4-rho*lna4;

f17=exp(lnc5)^(-1/gamma(5))*[1+phi*(exp(lnkp5)/exp(lnk5)-1)]-beta*exp(lncp5)^(-1/gamma(5))*[exp(lnthetap5)+exp(lnap5)*A*alpha*exp(lnkp5)^(alpha-1)*exp(lnlp5)^(1-alpha)];
f17b=b(5)*[exp(lnl5)^(alpha-1/eta(5))]-[exp(lnc5)^(-1/gamma(5))]*exp(lna5)*(1-alpha)*A*exp(lnk5)^alpha;
f18=exp(lntheta5)-[1+phi*[1+0.5*(exp(lnkp5)/exp(lnk5)-1)]*(exp(lnkp5)/exp(lnk5)-1)];
f19=tau(5)*exp(lnc5)^(-1/gamma(5))-exp(lnmu);
f20=lnap5-rho*lna5;

f21=exp(lnc6)^(-1/gamma(6))*[1+phi*(exp(lnkp6)/exp(lnk6)-1)]-beta*exp(lncp6)^(-1/gamma(6))*[exp(lnthetap6)+exp(lnap6)*A*alpha*exp(lnkp6)^(alpha-1)*exp(lnlp6)^(1-alpha)];
f21b=b(6)*[exp(lnl6)^(alpha-1/eta(6))]-[exp(lnc6)^(-1/gamma(6))]*exp(lna6)*(1-alpha)*A*exp(lnk6)^alpha;
f22=exp(lntheta6)-[1+phi*[1+0.5*(exp(lnkp6)/exp(lnk6)-1)]*(exp(lnkp6)/exp(lnk6)-1)];
f23=tau(6)*exp(lnc6)^(-1/gamma(6))-exp(lnmu);
f24=lnap6-rho*lna6;

f25=exp(lnc7)^(-1/gamma(7))*[1+phi*(exp(lnkp7)/exp(lnk7)-1)]-beta*exp(lncp7)^(-1/gamma(7))*[exp(lnthetap7)+exp(lnap7)*A*alpha*exp(lnkp7)^(alpha-1)*exp(lnlp7)^(1-alpha)];
f25b=b(7)*[exp(lnl7)^(alpha-1/eta(7))]-[exp(lnc7)^(-1/gamma(7))]*exp(lna7)*(1-alpha)*A*exp(lnk7)^alpha;
f26=exp(lntheta7)-[1+phi*[1+0.5*(exp(lnkp7)/exp(lnk7)-1)]*(exp(lnkp7)/exp(lnk7)-1)];
f27=tau(7)*exp(lnc7)^(-1/gamma(7))-exp(lnmu);
f28=lnap7-rho*lna7;

f29=exp(lnc8)^(-1/gamma(8))*[1+phi*(exp(lnkp8)/exp(lnk8)-1)]-beta*exp(lncp8)^(-1/gamma(8))*[exp(lnthetap8)+exp(lnap8)*A*alpha*exp(lnkp8)^(alpha-1)*exp(lnlp8)^(1-alpha)];
f29b=b(8)*[exp(lnl8)^(alpha-1/eta(8))]-[exp(lnc8)^(-1/gamma(8))]*exp(lna8)*(1-alpha)*A*exp(lnk8)^alpha;
f30=exp(lntheta8)-[1+phi*[1+0.5*(exp(lnkp8)/exp(lnk8)-1)]*(exp(lnkp8)/exp(lnk8)-1)];
f31=tau(8)*exp(lnc8)^(-1/gamma(8))-exp(lnmu);
f32=lnap8-rho*lna8;

f33=exp(lnc1)+exp(lnc2)+exp(lnc3)+exp(lnc4)+exp(lnc5)+exp(lnc6)+exp(lnc7)+exp(lnc8)-exp(lnC);
f34=exp(lnk1)+exp(lnk2)+exp(lnk3)+exp(lnk4)+exp(lnk5)+exp(lnk6)+exp(lnk7)+exp(lnk8)-exp(lnK);
f35=exp(lna1)*A*exp(lnk1)^alpha*exp(lnl1)^(1-alpha)+exp(lna2)*A*exp(lnk2)^alpha*exp(lnl2)^(1-alpha)+...
    exp(lna3)*A*exp(lnk3)^alpha*exp(lnl3)^(1-alpha)+exp(lna4)*A*exp(lnk4)^alpha*exp(lnl4)^(1-alpha)+...
    exp(lna5)*A*exp(lnk5)^alpha*exp(lnl5)^(1-alpha)+exp(lna6)*A*exp(lnk6)^alpha*exp(lnl6)^(1-alpha)+...
    exp(lna7)*A*exp(lnk7)^alpha*exp(lnl7)^(1-alpha)+exp(lna8)*A*exp(lnk8)^alpha*exp(lnl8)^(1-alpha)-exp(lnY);
f36=exp(lnC)+exp(lnKp)-exp(lnK)-exp(lnY)+...
    phi/2*[exp(lnk1)*(exp(lnkp1)/exp(lnk1)-1)^2+exp(lnk2)*(exp(lnkp2)/exp(lnk2)-1)^2+...
    exp(lnk3)*(exp(lnkp3)/exp(lnk3)-1)^2+exp(lnk4)*(exp(lnkp4)/exp(lnk4)-1)^2+...
    exp(lnk5)*(exp(lnkp5)/exp(lnk5)-1)^2+exp(lnk6)*(exp(lnkp6)/exp(lnk6)-1)^2+...
    exp(lnk7)*(exp(lnkp7)/exp(lnk7)-1)^2+exp(lnk8)*(exp(lnkp8)/exp(lnk8)-1)^2];

f=[f1;f1b;f2;f3;f4;f5;f5b;f6;f7;f8;f9;f9b;f10;f11;f12;f13;f13b;f14;f15;f16;f17;f17b;f18;f19;f20;f21;f21b;f22;f23;f24;f25;f25b;f26;f27;f28;f29;f29b;f30;f31;f32;f33;f34;f35;f36];

% 3) Define vectors of states and controls
xx=[lnk1 lnk2 lnk3 lnk4 lnk5 lnk6 lnk7 lnk8 lna1 lna2 lna3 lna4 lna5 lna6 lna7 lna8];
xxp=[lnkp1 lnkp2 lnkp3 lnkp4 lnkp5 lnkp6 lnkp7 lnkp8 lnap1 lnap2 lnap3 lnap4 lnap5 lnap6 lnap7 lnap8];

yy=[lnc1 lnc2 lnc3 lnc4 lnc5 lnc6 lnc7 lnc8 lntheta1 lntheta2 lntheta3 lntheta4 lntheta5 lntheta6 lntheta7 lntheta8 lnmu lnC lnK lnY lnl1 lnl2 lnl3 lnl4 lnl5 lnl6 lnl7 lnl8];
yyp=[lncp1 lncp2 lncp3 lncp4 lncp5 lncp6 lncp7 lncp8 lnthetap1 lnthetap2 lnthetap3 lnthetap4 lnthetap5 lnthetap6 lnthetap7 lnthetap8 lnmup lnCp lnKp lnYp lnlp1 lnlp2 lnlp3 lnlp4 lnlp5 lnlp6 lnlp7 lnlp8];

% 4) Compute analytical derivatives and evaluate them at the steady state
nx = length(xx);
ny = length(yy);

[fx,fxp,fy,fyp,fypyp,fypy,fypxp,fypx,fyyp,fyy,fyxp,fyx,fxpyp,fxpy,...
    fxpxp,fxpx,fxyp,fxy,fxxp,fxx]=anal_deriv(f,xx,yy,xxp,yyp);

assign_steadystate8;
num_eval;

% 5) Construct log-linear approximation using the Klein (2000) algorithm
[gx,hx] = gx_hx(nfy,nfx,nfyp,nfxp);



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 3. Construct initial guess for decision rule coefficients %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 1) Construct bounds for the state space
tr=1;  % scale parameter
kmin=1/1.05;
kmax=1.05;
amin=exp(-tr*2*sigma/sqrt(1-rho^2));
amax=exp(+tr*2*sigma/sqrt(1-rho^2));

 tr=0.8;  % scale parameter
 kmin=0.8;
 kmax=1.2;
 amin=exp(-tr*sigma/(1-rho));
 amax=exp(+tr*sigma/(1-rho));


% 2) Construct Smolyak Grid
nx = 2*N;
x_Smol = zeros(1,5);
for j = 1:5
  x_Smol(j) = -cos((j-1)*pi/4);
end

%Build Grid
xgrid = zeros(1+4*nx+2*nx*(nx-1),nx);
xgrid = x_Smol(3) * xgrid;
for j = 1:nx
  xgrid(1+4*(j-1)+1,j) = x_Smol(1);
  xgrid(1+4*(j-1)+2,j) = x_Smol(2);
  xgrid(1+4*(j-1)+3,j) = x_Smol(4);
  xgrid(1+4*(j-1)+4,j) = x_Smol(5);
end
ind = 1+4*nx+1;
points1 = [x_Smol(1);x_Smol(1);x_Smol(5);x_Smol(5)];
points2 = [x_Smol(1);x_Smol(5);x_Smol(1);x_Smol(5)];
for j = 1:(nx-1)
  for i = j+1: nx
    xgrid(ind:ind+3,j) = points1;
    xgrid(ind:ind+3,i) = points2;
    ind = ind+4;
  end
end


% Construct state vectors from grid
k1_vec=scalup(xgrid(:,1),kmin,kmax);
k2_vec=scalup(xgrid(:,2),kmin,kmax);
k3_vec=scalup(xgrid(:,3),kmin,kmax);
k4_vec=scalup(xgrid(:,4),kmin,kmax);
k5_vec=scalup(xgrid(:,5),kmin,kmax);
k6_vec=scalup(xgrid(:,6),kmin,kmax);
k7_vec=scalup(xgrid(:,7),kmin,kmax);
k8_vec=scalup(xgrid(:,8),kmin,kmax);

lna1_vec=scalup(xgrid(:,9),log(amin),log(amax));
lna2_vec=scalup(xgrid(:,10),log(amin),log(amax));
lna3_vec=scalup(xgrid(:,11),log(amin),log(amax));
lna4_vec=scalup(xgrid(:,12),log(amin),log(amax));
lna5_vec=scalup(xgrid(:,13),log(amin),log(amax));
lna6_vec=scalup(xgrid(:,14),log(amin),log(amax));
lna7_vec=scalup(xgrid(:,15),log(amin),log(amax));
lna8_vec=scalup(xgrid(:,16),log(amin),log(amax));

% Construct regressor matrix to be used for perturbation initialization
INIT_REG=[log(k1_vec)-log(kss) log(k2_vec)-log(kss) log(k3_vec)-log(kss) log(k4_vec)-log(kss) log(k5_vec)-log(kss) log(k6_vec)-log(kss) log(k7_vec)-log(kss) log(k8_vec)-log(kss)...
    lna1_vec-log(ass) lna2_vec-log(ass) lna3_vec-log(ass) lna4_vec-log(ass) lna5_vec-log(ass) lna6_vec-log(ass) lna7_vec-log(ass) lna8_vec-log(ass)];

% Construct next period's capital stocks and current consumption
lnkp1_vec=log(kss)+INIT_REG*hx(1,:)';
lnkp2_vec=log(kss)+INIT_REG*hx(2,:)';
lnkp3_vec=log(kss)+INIT_REG*hx(3,:)';
lnkp4_vec=log(kss)+INIT_REG*hx(4,:)';
lnkp5_vec=log(kss)+INIT_REG*hx(5,:)';
lnkp6_vec=log(kss)+INIT_REG*hx(6,:)';
lnkp7_vec=log(kss)+INIT_REG*hx(7,:)';
lnkp8_vec=log(kss)+INIT_REG*hx(8,:)';
kp1_vec = exp(lnkp1_vec);
kp2_vec = exp(lnkp2_vec);
kp3_vec = exp(lnkp3_vec);
kp4_vec = exp(lnkp4_vec);
kp5_vec = exp(lnkp5_vec);
kp6_vec = exp(lnkp6_vec);
kp7_vec = exp(lnkp7_vec);
kp8_vec = exp(lnkp8_vec);

lnc1_vec=lncss+INIT_REG*gx(1,:)';
lnc2_vec=lncss+INIT_REG*gx(2,:)';
lnc3_vec=lncss+INIT_REG*gx(3,:)';
lnc4_vec=lncss+INIT_REG*gx(4,:)';
lnc5_vec=lncss+INIT_REG*gx(5,:)';
lnc6_vec=lncss+INIT_REG*gx(6,:)';
lnc7_vec=lncss+INIT_REG*gx(7,:)';
lnc8_vec=lncss+INIT_REG*gx(8,:)';
c1_vec = exp(lnc1_vec);
c2_vec = exp(lnc2_vec);
c3_vec = exp(lnc3_vec);
c4_vec = exp(lnc4_vec);
c5_vec = exp(lnc5_vec);
c6_vec = exp(lnc6_vec);
c7_vec = exp(lnc7_vec);
c8_vec = exp(lnc8_vec);

OUT_VEC=[kp1_vec c1_vec kp2_vec c2_vec kp3_vec c3_vec kp4_vec c4_vec ...
        kp5_vec c5_vec kp6_vec c6_vec kp7_vec c7_vec kp8_vec c8_vec];
for j = 1:(1+4*nx+2*nx*(nx-1))
  for i = 1:N
  OUT_VEC1(N*(j-1)+i,:) = OUT_VEC(j,2*i-1:2*i);
  end
end

% Output linearization to .csv file
dlmwrite('initguess8.csv',OUT_VEC1,'delimiter',',','precision',12)



