% Program that solves Greenwood et. all (1988) model with Capacity utilization
% via quadratic aproximation
%
% by Boubacar Drame,
%
% Nottingham, April 05, 2010 
% NB: For more detailed comments on the codes below please refer to the
% appendix section of the paper submitted along with the present.

clear all;
clc;

%parameter values
beta = .96;
alpha = .29;
gamma= 2;
theta = .6;
omega= 1.42;
rho=.51;
sigmaeps = .045;

%Steady state equations.

xss=0;
xssprime=xss;
hss=((1-beta)/beta*(omega/(omega-1)))^(1/omega)
kss=(1-alpha)^(1/theta)*alpha^((alpha+theta)/((1-alpha)*theta))*hss^((alpha*(1+theta)-omega*(alpha+theta))/(theta*(1-alpha)))
lss=(1-alpha)^(1/(alpha+theta))*kss^(alpha/(alpha+theta))*hss^(alpha/(alpha+theta))
css=((kss*hss)^alpha)*lss^(1-alpha)-(1/omega)*(hss^omega)*kss
kssprime=kss;

h=log(((1-beta)/beta*(omega/(omega-1)))^(1/omega))

sss = [kss; xss; lss; hss; kssprime]; %vector of steady state values

% This is our Objective function for the Quadratic Approximation

fun = @utilfun;

% QA for the utility function
U = utilfun(sss);
J = gradp(fun,sss);     % Jacobian evaluated in ss
H = hessp(fun,sss);     % Hessian evaluated in ss
r = sss;

% Where we can form Qup matrix
Qu11 = U-r'*J+0.5*r'*H*r;
Qu12 = 0.5*(J-H*r)';
Qu21 = 0.5*(J-H*r);
Qu22 = 0.5*H;
Qu = [Qu11 Qu12; Qu21 Qu22];  % 6x6 Qu matrix -> r' Qu r

V = -.01*eye(3); % A good first guess is the negative definite 3*3 identity matrix


%/*Order of elements associated with W*/

%/*	1	k	x	l   h   kprime   1	 k'	x'  */

%/*	1	2	3	4	5	6	     7   8  9   */


% We now form the H matrix to impose restrictions and thus having an
% unconstrained optimisation problem

H = zeros(9,6);

H(1:6,1:6) = eye(6);
H(7,1)= 1;   % 1 = 1
H(8,6)=1;     % kprime = k' 
H(9,3)=rho;  % Restriction: x' = rho * x


tdif = 10; 	%/*Convergence criterion*
% This is set arbitrarilly large in order to enter the iterative loop

while tdif>.0001;

W = zeros(9,9);	%/*Formation of the large matrix for the RHS of the Bellman equation
W(1:6,1:6) = Qu;
W(7:9,7:9)=beta*V;

W6 = H'*W*H;		%/We now impose the constraints

B = [-(inv(W6(4:6,4:6))*W6(4:6,1:3))];

%The H1 matrix is composed of the identity matrix (on top)and B at the
%bottom

H1= [eye(3); B];

W3=H1'*W6*H1;	%/*Derivation of the update for the value function through imposing the decision rule of investment

dif = abs(W3-V);

tdif = sum(sum(dif));


	V=W3;

% tdif 		%/*Help keep track of the convergence process

end;

%/*The respective decision rules are 

lc=B(1,1);
lk=B(1,2);
lx=B(1,3);

hc=B(2,1);
hk=B(2,2);
hx=B(2,3);

kprimec=B(3,1);
kprimek=B(3,2);
kprimex=B(3,3);


%/******************************************************************************/
%/ Performing Monte Carlo simulation 
%/******************************************************************************/

m=100;		%/Number of experiments*/

n = 10000;	%/Length of series for each experiment*/


%/*Initialize the Monte-Carlo Simulation vectors*/

stdy = zeros(m,1);
stdcons = zeros(m,1);
stdinvestm = zeros(m,1);
correlycm = zeros(m,1);
correlyim = zeros(m,1);
correlym = zeros(m,1);
correlcm = zeros(m,1);
correlim = zeros(m,1);
csym = zeros(m,1);



for i=1:m; % Start loop for m


k0 = kss;  %/*starting value for capital*/

k=zeros(n,1);
kdash=zeros(n+1,1);
x=zeros(n+1,1);
eps=zeros(n+1,1);
y=zeros(n,1);
c=zeros(n,1);
investm=zeros(n,1);
l=zeros(n,1);
h=zeros(n,1);

k(1)=k0;
eps(1)= sigmaeps*randn(1);
x(1) = eps(1);


for j=1:n;  %Start loop for n

eps(j+1)= sigmaeps*randn(1,1);
x(j+1)=x(j)*rho+eps(j+1);
kprime(j)= kprimec + kprimek*k(j) + kprimex*x(j);
l(j)= lc + lk*k(j) + lx*x(j);
h(j) = hc + hk*k(j) + hx*x(j);
y(j) = ((h(j)*k(j))^alpha)*(l(j)^(1-alpha));
investm(j) = exp(-x(j))*kprime(j) - exp(-x(j))*k(j)*(1-((h(j)^omega)/omega));
c(j)=y(j)-investm(j);
k(j+1)=kprime(j);

end; % End loop for n


lnc = log(c);
lny = log(y);
lnh = log(h);
lnl = log(l);
lninvestm = log(investm);

%reducing the n+1 vector to n elements

lnkhelp = log(k);
lnk=lnkhelp(1:n,1);

%For Productivity  we have output divided by total hours, using logs.

lnp = lny - lnl;


%Calculate correlation between each output and each variable of interest

correlyc = corrcoef(lnc, lny);
correlyi = corrcoef(lninvestm, lny);
correlyy = corrcoef(lny, lny);
correlyh = corrcoef(lnh, lny);
correlyl = corrcoef(lnl, lny);
correlyk = corrcoef(lnk, lny);
correlyp = corrcoef(lnp, lny);

correlycm(i) = correlyc(2,1);
correlyim(i) = correlyi(2,1);
correlyym(i) = correlyy(2,1);
correlyhm(i) = correlyh(2,1);
correlylm(i) = correlyl(2,1);
correlykm(i) = correlyk(2,1);
correlypm(i) = correlyp(2,1);

%Calculate standatd deviation for each variable of interest

stdy(i) = 100*std(lny);
stdcons(i) = 100*std(lnc);
stdinvestm(i) = 100*std(lninvestm);
stdutil(i) = 100*std(lnh);
stdlabour(i) = 100*std(lnl);
stdcapital(i) = 100*std(lnk);
stdproductivity(i) = 100*std(lnp);

%Calculate autocorrelation for each variable of interest


correly=corrcoef(lny(1:(length(lny)-1)), lny(2:(length(lny))));
correlc=corrcoef(lnc(1:(length(lnc)-1)), lnc(2:(length(lnc))));
correli=corrcoef(lninvestm(1:(length(lninvestm)-1)), lninvestm(2:(length(lninvestm))));
correlh=corrcoef(lnh(1:(length(lnh)-1)), lnh(2:(length(lnh))));
correll=corrcoef(lnl(1:(length(lnl)-1)), lnl(2:(length(lnl))));
correlk=corrcoef(lnk(1:(length(lnl)-1)), lnk(2:(length(lnl))));
correlp=corrcoef(lnp(1:(length(lnl)-1)), lnp(2:(length(lnl))));

correlym(i) = correly(2,1);
correlcm(i) = correlc(2,1);
correlim(i) = correli(2,1);
correlhm(i) = correlh(2,1);
correllm(i) = correll(2,1);
correlkm(i) = correlk(2,1);
correlpm(i) = correlp(2,1);


end; %loop for m


% Report Monte-Carlo results
char 'std(Output), Paper = 3.5:'
mean(stdy)

char 'std(Consumption), Paper = 2.2:'
mean(stdcons)

char 'std(Investment), Paper = 11.6:'
mean(stdinvestm)

char 'std(Utilisation Rate), Paper = 6.0:'
mean(stdutil)

char 'std(Hours), Paper = 2.2:'
mean(stdlabour)

char 'std(Capital Stock), Paper = 5.6:'
mean(stdcapital)

char 'std(Productivity), Paper = 1.3:'
mean(stdproductivity)


char 'Corr(Output,Output):, Paper = 1.0'
mean(correlyym)

char 'Corr(Consumption,Output), Paper = 0.79:'
mean(correlycm)

char 'Corr(Investment,Output):, Paper = 0.90'
mean(correlyim)

char 'Corr(Utilisation Rate,Output), Paper = 0.61:'
mean(correlyhm)

char 'Corr(Hours,Output), Paper = 1.00:'
mean(correlylm)

char 'Corr(Capital Stock,Output), Paper = 0.52:'
mean(correlykm)

char 'Corr(Productivity,Output), Paper = 1.00:'
mean(correlypm)



char 'autocorrelation Output, Paper = 0.66:'
mean(correlym)

char 'autocorrelation Consumption, Paper = 0.94:'
mean(correlcm)

char 'autocorrelation Investment, Paper = 0.50:'
mean(correlim)

char 'autocorrelation Utilisation Rate, Paper = 0.52:'
mean(correlhm)

char 'autocorrelation Hours, Paper = 0.66:'
mean(correllm)

char 'autocorrelation Capital Stock, Paper = 0.99:'
mean(correlkm)

char 'autocorrelation Productivity, Paper = 0.66:'
mean(correlpm)