%% 14.461 Problem Set 1 Fall 2009
% Question 1

clear all
close all
clc

%% Parameters

b = 0.99; d = 0.025; e =1; sg = 1; a = 0.36; ro = 0.95; A = 1;

p1 = ((1/b -1 + d) / (A*a)) ^ (1/ (1-a));
p2 = A * (p1) ^ (-a) - d * (p1)^(-1);

%% Steady State

N = ((1-a)*A * (p2) ^(-sg) * (p1) ^ (-a) ) ^ (1 / (e+sg));
C = p2 * N;
K = N / p1;
Y = A * K ^ a * N ^ (1-a);

%% Define more parameters

p3 = (1/b * K / C  + (1-a) * Y /C * (e+a)^ (-1) * a  )/ (1 + (1-a) * Y / C * (e+a)^ (-1) * sg   ) ;
p4 = (K / C) / ( 1 + (1-a) * Y / C * (e+a)^ (-1) * sg );
p5 = ( (1-a) * Y / C * (e+a)^ (-1) + Y / C ) / ( 1 + (1-a) * Y / C * (e+a)^ (-1) * sg );
p6 = (K / C) / ( 1 + (1-a) * Y / C * (e+a)^ (-1) * sg );

%% Main quadratic equation

aa = p4 * (1 + (1 - (1 - d )*b )*(1-a) * sg / (e+a) );
bb = -p3 - sg*p4  +  (1 - (1 - d )*b )*(1-a) * ( (a - sg * p3) / (e+a) - 1)  ;
cc = sg * p3;

vkk = (-bb - sqrt (bb^2 - 4 * aa * cc)) / (2 * aa);
vck = p3 - p4 * vkk;
vnk = a / (a + e) - sg / (a + e) * p3 + (sg / (a + e)) * p4 * vkk;

vka =  - ( (1 - (1 - d )*b )*ro + (1 - (1 - d )*b )*(1-a)*ro / (e+a) *(1 - sg*p5 ) + sg * (1-ro) * p5 ) / ( (1 - (1 - d )*b )*(1-a) * ( ro * p6 * sg /(e+a) - 1 ) - sg * (1 - ro) * p6 - vck + vnk * (1 - (1 - d )*b )*(1-a) );
vca = p5 - p6 * vka;
vna = 1 / (e+a)  - sg / (e+a) * p5 + sg / (e+a )*p6 *vka;

%% IR functions

eps(1) = 1; eps(2:50) = zeros(49,1);

aaa(1) = 0;
kkk(1) = 0;
ccc(1) = 0;
nnn(1) = 0;
yyy(1) = 0;
iii(1) = 0

for i=1:50
    aaa(i+1) = ro * aaa(i) + eps(i);
    kkk(i+1) = vkk * kkk(i) + vka * aaa(i+1);
    nnn(i+1) = vnk * kkk(i) + vna * aaa(i+1);
    ccc(i+1) = vck * kkk(i) + vca * aaa(i+1);
    yyy(i+1) = aaa(i+1) + a * kkk(i) + (1-a) * nnn(i+1);
    
end



plot(ccc)