%%%%% Problem set 1 
%%%%% Declare parameter values

alpha = 0.33;
delta = 0.04;
z     = 1;
sigma = 2; 
beta  = 0.97;

parameters1 = [z,alpha,delta,beta];

%%%%% Declare some symbols

syms kk cc zz aa dd bb 

[cbar,kbar] = solve('cc+kk-zz*(kk^aa)-(1-dd)*kk=0',...
                    '1-bb*1*(1+aa*zz*(kk^(aa-1))-dd)=0',cc,kk); 
                % Note: The second equation can only be solved for kk since 
                % cc does not appear in it. 

kbar1 = double(subs(kbar,{zz,aa,dd,bb},parameters1));           
cbar1 = double(subs(cbar,{zz,aa,dd,bb},parameters1));

ybar1 = 1*(kbar1^alpha); %% output
ibar1 = delta*kbar1;     %% investment
wbar1 = (1-alpha)*ybar1; %% wage rate
rbar1 = alpha*ybar1/kbar1; %% rental rate
                       
%%%%% KNOWN COEFFICIENTS

A = 1;
B = -1/beta;
C = cbar1/kbar1;
F = (1-alpha)*((1-beta+beta*delta)/sigma);
G = 0;
J = 1;
K = -1;

%%%%% SOLVE QUADRATIC IN P

quadratic_coefficients = [-J*inv(C)*A,F-J*inv(C)*B-K*inv(C)*A,G-K*inv(C)*B];

Ps = roots(quadratic_coefficients) %% "roots" solves polynomials

P  = min(Ps) %% check that this is real and less than one in modulus.

R = -inv(C)*(A*P+B) %% back out associated value of R

%%%%% QUESTION 3 and QUESTION5a
T  = 100;
k0 = -0.10;

ks = zeros(T,1); %% vectors to store answers
cs = zeros(T,1);
ys = zeros(T,1);
is = zeros(T,1);
ws = zeros(T,1);
rs = zeros(T,1);

ts = [0;cumsum(ones(T-1,1))];

ks(1)=k0;

for t = 1:T-1,
    ks(t+1) = P*ks(t);
    cs(t)   = R*ks(t);
    ys(t)   = alpha*ks(t);
    is(t)   = inv(delta)*(ks(t+1)-(1-delta)*ks(t));
    ws(t)   = alpha*ks(t);
    rs(t)   = -(1-alpha)*ks(t);
end

figure(3)
subplot(2,1,1)
plot(ts,ks,ts,cs,ts,ys,ts,is,ts,zeros(T,1),'k--')
title('Recovery from a war')
ylabel('Proportional deviation')
legend('capital','consumption','output','investment',4)
subplot(2,1,2)
plot(ts,ws,ts,rs,ts,zeros(T,1),'k--')
ylabel('Proportional deviation')
legend('wages','rental rate')

%%%%% QUESTION 4 and QUESTION5b

%%    First compute new steady state

parameters2 = [2,alpha,delta,beta];

kbar2 = double(subs(kbar,{zz,aa,dd,bb},parameters2));           
cbar2 = double(subs(cbar,{zz,aa,dd,bb},parameters2));
ybar2 = 2*(kbar2^alpha); %% output
ibar2 = delta*kbar2;     %% investment
wbar2 = (1-alpha)*ybar2; %% wage rate
rbar2 = alpha*ybar2/kbar2; %% rental rate

%%    Now compute implied deviation relative to new steady state

k0implied = log(kbar1/kbar2);

T  = 100;

ks = zeros(T,1); %% vectors to store answers
cs = zeros(T,1);
ys = zeros(T,1);
is = zeros(T,1);
ws = zeros(T,1);
rs = zeros(T,1);

ks(1)=k0implied;

%%  Compute time paths relative to new steady state

for t = 1:T-1,
    ks(t+1) = P*ks(t);
    cs(t)   = R*ks(t);
    ys(t)   = alpha*ks(t);
    is(t)   = inv(delta)*(ks(t+1)-(1-delta)*ks(t));
    ws(t)   = alpha*ks(t);
    rs(t)   = -(1-alpha)*ks(t);
     
end

%%   Now measure everything relative to old steady state values

ks = ks+log(kbar2/kbar1);
cs = cs+log(cbar2/cbar1);
ys = ys+log(ybar2/ybar1);
is = is+log(ibar2/ibar1);
ws = ws+log(wbar2/wbar1);
rs = rs+log(rbar2/rbar1);

%%   Now include a few t<0 for the pictures

ks = [0;0;0;0;ks];
cs = [0;0;0;0;cs];
ys = [0;0;0;0;ys];
is = [0;0;0;0;is];
ws = [0;0;0;0;ws];
rs = [0;0;0;0;rs];
ts = [-4;-3;-2;-1;0;cumsum(ones(T-1,1))];

figure(4)
subplot(2,1,1)
plot(ts,ks,ts,cs,ts,ys,ts,is,ts,zeros(length(ts),1),'k--')
title('Technological breakthrough')
ylabel('Proportional deviation')
legend('capital','consumption','output','investment',4)
subplot(2,1,2)
plot(ts,ws,ts,rs,ts,zeros(length(ts),1),'k--')
ylabel('Proportional deviation')
legend('wages','rental rate')