% Main file to generate initial guess for solution to model A5 
% (asymmetric CRRA utility function 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
global kmin kmax amin amax
global xgrid
global k1_mat k2_mat k3_mat k4_mat k5_mat k6_mat k7_mat k8_mat k9_mat k10_mat K_mat
global a1_mat a2_mat a3_mat a4_mat a5_mat a6_mat a7_mat a8_mat a9_mat a10_mat
global ap1_mat ap2_mat ap3_mat ap4_mat ap5_mat ap6_mat ap7_mat ap8_mat ap9_mat ap10_mat
global y1_mat y2_mat y3_mat y4_mat y5_mat y6_mat y7_mat y8_mat y9_mat y10_mat Y_mat

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

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

% 1) Calibration
N=10;
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];

% 2) Compute the steady state and the Negishi weights
kss=1;
Kss=N*kss;
css=A*kss^alpha;
Css=N*css;
lnass=0;
ass=1;
thetass=1;
Yss=N*A*kss^alpha;


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


% 4) Compute the steady state in logs
lnkss=log(kss);
lnKss=log(Kss);
lncss=log(css);
lnCss=log(Css);
lnthetass=log(thetass);
lnYss=log(Yss);
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
syms c2 lna2 k2 theta2 cp2 lnap2 kp2 thetap2
syms c3 lna3 k3 theta3 cp3 lnap3 kp3 thetap3
syms c4 lna4 k4 theta4 cp4 lnap4 kp4 thetap4
syms c5 lna5 k5 theta5 cp5 lnap5 kp5 thetap5
syms c6 lna6 k6 theta6 cp6 lnap6 kp6 thetap6
syms c7 lna7 k7 theta7 cp7 lnap7 kp7 thetap7
syms c8 lna8 k8 theta8 cp8 lnap8 kp8 thetap8
syms c9 lna9 k9 theta9 cp9 lnap9 kp9 thetap9
syms c10 lna10 k10 theta10 cp10 lnap10 kp10 thetap10
syms C Cp K Kp Y Yp CAC CACp mu mup

syms lnc1 lna1 lnk1 lntheta1 lncp1 lnap1 lnkp1 lnthetap1
syms lnc2 lna2 lnk2 lntheta2 lncp2 lnap2 lnkp2 lnthetap2
syms lnc3 lna3 lnk3 lntheta3 lncp3 lnap3 lnkp3 lnthetap3
syms lnc4 lna4 lnk4 lntheta4 lncp4 lnap4 lnkp4 lnthetap4
syms lnc5 lna5 lnk5 lntheta5 lncp5 lnap5 lnkp5 lnthetap5
syms lnc6 lna6 lnk6 lntheta6 lncp6 lnap6 lnkp6 lnthetap6
syms lnc7 lna7 lnk7 lntheta7 lncp7 lnap7 lnkp7 lnthetap7
syms lnc8 lna8 lnk8 lntheta8 lncp8 lnap8 lnkp8 lnthetap8
syms lnc9 lna9 lnk9 lntheta9 lncp9 lnap9 lnkp9 lnthetap9
syms lnc10 lna10 lnk10 lntheta10 lncp10 lnap10 lnkp10 lnthetap10
syms lnC lnCp lnK lnKp lnY lnYp CAC CACp lnmu lnmup

% 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)];
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)];
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)];
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)];
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)];
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)];
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)];
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)];
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(lnc9)^(-1/gamma(9))*[1+phi*(exp(lnkp9)/exp(lnk9)-1)]-beta*exp(lncp9)^(-1/gamma(9))*[exp(lnthetap9)+exp(lnap9)*A*alpha*exp(lnkp9)^(alpha-1)];
f34=exp(lntheta9)-[1+phi*[1+0.5*(exp(lnkp9)/exp(lnk9)-1)]*(exp(lnkp9)/exp(lnk9)-1)];
f35=tau(9)*exp(lnc9)^(-1/gamma(9))-exp(lnmu);
f36=lnap9-rho*lna9;

f37=exp(lnc10)^(-1/gamma(10))*[1+phi*(exp(lnkp10)/exp(lnk10)-1)]-beta*exp(lncp10)^(-1/gamma(10))*[exp(lnthetap10)+exp(lnap10)*A*alpha*exp(lnkp10)^(alpha-1)];
f38=exp(lntheta10)-[1+phi*[1+0.5*(exp(lnkp10)/exp(lnk10)-1)]*(exp(lnkp10)/exp(lnk10)-1)];
f39=tau(10)*exp(lnc10)^(-1/gamma(10))-exp(lnmu);
f40=lnap10-rho*lna10;

f41=exp(lnc1)+exp(lnc2)+exp(lnc3)+exp(lnc4)+exp(lnc5)+exp(lnc6)+exp(lnc7)+exp(lnc8)+exp(lnc9)+exp(lnc10)-exp(lnC);
f42=exp(lnk1)+exp(lnk2)+exp(lnk3)+exp(lnk4)+exp(lnk5)+exp(lnk6)+exp(lnk7)+exp(lnk8)+exp(lnk9)+exp(lnk10)-exp(lnK);
f43=exp(lna1)*A*exp(lnk1)^alpha+exp(lna2)*A*exp(lnk2)^alpha+exp(lna3)*A*exp(lnk3)^alpha+exp(lna4)*A*exp(lnk4)^alpha+...
	exp(lna5)*A*exp(lnk5)^alpha+exp(lna6)*A*exp(lnk6)^alpha+exp(lna7)*A*exp(lnk7)^alpha+exp(lna8)*A*exp(lnk8)^alpha+exp(lna9)*A*exp(lnk9)^alpha+exp(lna10)*A*exp(lnk10)^alpha-exp(lnY);
f44=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+exp(lnk9)*(exp(lnkp9)/exp(lnk9)-1)^2+exp(lnk10)*(exp(lnkp10)/exp(lnk10)-1)^2];

f=[f1;f2;f3;f4;f5;f6;f7;f8;f9;f10;f11;f12;f13;f14;f15;f16;f17;f18;f19;f20;f21;f22;f23;f24;f25;f26;f27;f28;f29;f30;f31;f32;f33;f34;f35;f36;f37;f38;f39;f40;f41;f42;f43;f44];

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

yy=[lnc1 lnc2 lnc3 lnc4 lnc5 lnc6 lnc7 lnc8 lnc9 lnc10 lntheta1 lntheta2 lntheta3 lntheta4 lntheta5 lntheta6 lntheta7 lntheta8 lntheta9 lntheta10 lnmu lnC lnK lnY];
yyp=[lncp1 lncp2 lncp3 lncp4 lncp5 lncp6 lncp7 lncp8 lncp9 lncp10 lnthetap1 lnthetap2 lnthetap3 lnthetap4 lnthetap5 lnthetap6 lnthetap7 lnthetap8 lnthetap9 lnthetap10 lnmup lnCp lnKp lnYp];

% 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_steadystate10;
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);
k9_vec=scalup(xgrid(:,9),kmin,kmax);
k10_vec=scalup(xgrid(:,10),kmin,kmax);

lna1_vec=scalup(xgrid(:,11),log(amin),log(amax));
lna2_vec=scalup(xgrid(:,12),log(amin),log(amax));
lna3_vec=scalup(xgrid(:,13),log(amin),log(amax));
lna4_vec=scalup(xgrid(:,14),log(amin),log(amax));
lna5_vec=scalup(xgrid(:,15),log(amin),log(amax));
lna6_vec=scalup(xgrid(:,16),log(amin),log(amax));
lna7_vec=scalup(xgrid(:,17),log(amin),log(amax));
lna8_vec=scalup(xgrid(:,18),log(amin),log(amax));
lna9_vec=scalup(xgrid(:,19),log(amin),log(amax));
lna10_vec=scalup(xgrid(:,20),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) log(k9_vec)-log(kss) log(k10_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) lna9_vec-log(ass) lna10_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,:)';
lnkp9_vec=log(kss)+INIT_REG*hx(9,:)';
lnkp10_vec=log(kss)+INIT_REG*hx(10,:)';
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);
kp9_vec = exp(lnkp9_vec);
kp10_vec = exp(lnkp10_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,:)';
lnc9_vec=lncss+INIT_REG*gx(9,:)';
lnc10_vec=lncss+INIT_REG*gx(10,:)';
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);
c9_vec = exp(lnc9_vec);
c10_vec = exp(lnc10_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 kp9_vec c9_vec kp10_vec c10_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('initguess9.csv',OUT_VEC1,'delimiter',',','precision',12)
