%
%     peapro1.m MATLAB program to solve the neoclassical growth model
%
%     This program needs the following functions:
%
%     1. aexp.m
%     2. chebnode.m
%     3. chebpol.m
%     4. confun.m
%     5. gridfun.m
%     6. hernodes.m
%     7. makepoly.m
%     8. numi.m
%     9. savfun.m
%     10. scaldown.m
%     11. scalup.m
%     12. nsecant.m  
%
%     These files are available in the executable zip file pealib.exe
%
%
%              - Projection method & nonlinear equation solver to find
%                coefficients in approximating function
%
%      NOVEMBER 9 1998
%
%      Written by Christian Haefke, University of California, San Diego
%      Support by Wouter den Haan's NSF grant is gratefully acknowledged.
%
%      Documentation of the program is given at the end of the program
%
%
%
%============================================================================

clear all;
global alpha dfactor nu delta rho sigma po_k po_t tmin tmax kmin kmax;
global maxiter psitol;

!rm growth.out
diary growth.out

% 1. Initialize the parameters
% ============================

% 1.1 Initialize model parameters
% -------------------------------

alpha   = 0.33;         % Capital share 
dfactor = .95;          % Time discount factor
nu      = 2.;            % Risk aversion parameter
delta   = .95;          % Depreciation rate

sigma   = .010;         % Standard Deviation for log noise in technology shock
rho     = 0.95;         % Persistence of log technology shock

simyes  = 1;            % Simulate the model after solution
simtime = 1000;          % Length of simulation

graphs  = 1;            % Plot graphs for simulation
accur   = 1;            % Perform accuracy check

% 1.2 Initialize algorithm parameters
% -----------------------------------

po_k    =  3;           % Order of Polynomial for k
po_t    =  3;           % Order of Polynomial for theta

knod    =  po_k+1;           % Gridpoints for k
tnod    =  po_t+1;           % Gridpoints for theta

kmin    = 0.8;         % Lowest  gridpoint for capital 
			% as a percentage of steady state
kmax    = 1.2;         % Highest gridpoint for capital 
			% as a percentage of steady state

tmin    =-3.5;          % Lowest  gridpoint for log theta 
			% in standard deviations away from mean
tmax    = 3.5;          % Highest gridpoint for log theta 
			% in standard deviations away from mean

psi0file= 'psi0.in';   % File containig initial values for polynomial

maxiter = 2000;         % Maximum Iterations to find
			% parameters of polynomial
psitol  = 1e-5;         % Convergence criterion
			

% 1.3 Compute the Steady State
% ----------------------------

ks      = ( (1-dfactor+delta*dfactor) / (alpha*dfactor) )^(1/(alpha-1));
kmin    = log(ks * kmin);
kmax    = log(ks * kmax);

tsig    = sigma/sqrt(1-rho*rho);
tmin    = tmin*tsig;
tmax    = tmax*tsig;


% ***************************************************************************
% *                        M a i n   P r o g r a m                          *
% ***************************************************************************


% 2 Some Preparatory Work
% =======================

% 2.1 Create the grid  
% -------------------
% All gridpoints are stored in a matrix xgrid with dimensions (knod*tnod) x 2.
% Note that the gridpoints are just the standard chebyshev nodes which are in the 
% interval [-1,1]. The function scalup will tell you what the value of the state variable
% is for any point in [-1,1].

kg      = chebnode(knod);
tg      = chebnode(tnod);

for ki  = 1:knod;
	for ti  = 1:tnod;
		xgrid((ki-1)*tnod+ti,:) = [kg(ki) tg(ti)];
	end;
end;

[dimgrid nvars] = size(xgrid);


% 2.3 Initialize the coefficients of the polynomial  (Step I)
% -----------------------------------------------------------

% 2.3.1 Compute the size
% ----------------------

psisize = (po_k+1)*(po_t+1);
psi     = zeros(psisize,1);

% 2.3.2. Load the file
% --------------------

load(psi0file);

% 2.3.3. Check the size
% ---------------------

if any(size(psi0) ~= size(psi));
	msg = sprintf(' must contain a column vector with%2.0f rows and 1 column \n not%2.0f rows and%2.0f columns!',psisize,size(psi0));
	error([psi0file msg]);
end;

% 2.3.4. Initialize psi
% ---------------------

psi     = psi0;

% ----------------------------------------------------------------------

% 3 Call the equation solver
% ==========================

[psi,F] = nsecant('gridfun',psi0,0,xgrid);

%
%       3.3 Print Parameters
%       --------------------

sum1    = sprintf('alpha \tdfactor\t nu    \t delta \t');
sum2    = sprintf('%4.2f \t%4.2f \t %4.2f \t %4.2f \t',alpha,dfactor,nu,delta);
sum3    = sprintf('rho \t sigma    \n' );
sum4    = sprintf('%4.2f \t %4.2f \n',rho,sigma);
disp([sum1 sum3 sum2 sum4]);

diary psi0.out
psiout  = sprintf('%8.7f \n',psi);
disp(psiout);
diary off

% ***************************************************************************
% *                          S I M U L A T I O N                            *
% ***************************************************************************

if simyes == 1;

% 4.0 Initialize matrices
% -----------------------

lncs      = zeros(simtime,1);
lnks      = zeros(simtime+1,1);
lnts      = zeros(simtime+1,1);
epsi      = randn(simtime+1,1) * sigma;

% 4.1 Start at the Steady State
% -----------------------

lnks(1)   = log(( (1-dfactor+delta*dfactor) / (alpha*dfactor) )^(1/(alpha-1)) );
lnts(1)   = 0;
ss      = sprintf('Steady State at capital %d ',lnks(1));
disp(ss);

% 4.1 Run the simulation
% ----------------------

for ti  = 1:simtime
	lncs(ti)  = confun(lnks(ti),lnts(ti),psi);
	lnks(ti+1)= savfun(lnks(ti),lnts(ti),psi);
	lnts(ti+1)= rho * lnts(ti) + epsi(ti+1);
end;

% 4.2 Evaluate the simulation
% ---------------------------
ev1     = sprintf('\n Over the simulation the means and variances were \n');
ev2     = sprintf('\t Capital \t Consumption \t Shock \n');
ev3     = sprintf('Mean \t %f \t %f \t %f \n',mean(lnks),mean(lncs),mean(lnts));
ev4     = sprintf('Std.Dev\t %f \t %f \t %f \n',cov(lnks)^.5,cov(lncs)^.5,cov(lnts)^.5);
ev5     = sprintf('Min \t %f \t %f \t %f \n',min(lnks),min(lncs),min(lnts));
ev6     = sprintf('Max \t %f \t %f \t %f \n',max(lnks),max(lncs),max(lnts));

disp([ev1 ev2 ev3 ev4 ev5 ev6]);

if graphs       == 1;
figure(1)
	h1=plot([lncs(1:simtime) lnks(1:simtime)]); 
	legend('Consumption','Capital Stock',0);
	title('Log of Consumption and Capital Stock Over the Simulation');
	xlabel('Time');
	ylabel('Consumption, Capital');
print -noui -deps2 -r300 simulation.eps;
end;
end;

% 4.3 Check the Accuracy
% ----------------------

if accur       == 1;
tobs    = ceil(tnod/2);
obs     = 50;
lnk1    = (kmin:(kmax-kmin)/(obs-1):kmax)';
lnt1    = scalup(tg(tobs)*(ones(obs,1)),tmin,tmax);
lnk2    = savfun(lnk1,lnt1,psi);
for ai  = 1:obs;
	actk(ai,1)= numi('aexp',lnk2(ai),lnt1(ai),psi);
end;
k1ls    = scaldown(lnk1,kmin,kmax);
t1ls    = scaldown(lnt1,tmin,tmax);
apoly   = makepoly([po_k po_t],[k1ls t1ls]);
appk    = exp(apoly*psi);

figure(2)
	h2=plot(lnk1,[actk appk]);
	set(h2,{'LineStyle'},{'-';'--'});
	set(h2,{'Color'},{'k';'r'});
	titel=sprintf('Actual and Approximate Expectations for ln \\theta_{t}=%1.2f',lnt1(1));
	title(titel);
	leg1=('Actual Expectation');
	leg2=('Approximate Expectation');
	legend(leg1,leg2);
	xlabel('Log of Initial Capital Stock, ln k_{t}');
	ylabel('\beta E(c_{t+1}^{-\gamma}(\alpha \theta_{t+1}k_{t+1}^{\alpha-1}+1 - \delta)|\theta_t)');
print -noui -deps2 -r300 accuracyk.eps;

kobs    = ceil(knod/2);
lnt1    = (tmin:(tmax-tmin)/(obs-1):tmax)';
lnk1    = scalup(kg(kobs)*ones(obs,1),kmin,kmax);
lnk2    = savfun(lnk1,lnt1,psi);
for ai  = 1:obs;
	actt(ai,1)= numi('aexp',lnk2(ai),lnt1(ai),psi);
end;
k1ls    = scaldown(lnk1,kmin,kmax);
t1ls    = scaldown(lnt1,tmin,tmax);
apoly   = makepoly([po_k po_t],[k1ls t1ls]);
appt    = exp(apoly*psi);

figure(3)
	h3=plot(lnt1,[actt appt]);
	set(h3,{'LineStyle'},{'-';'--'});
	set(h3,{'Color'},{'k';'r'});
	titel=sprintf('Actual and Approximate Expectations for ln k_t=%1.2f',lnk1(1));
	title(titel);
	leg1=sprintf('Actual Expectation');
	leg2=sprintf('Approximate Expectation');
	legend(leg1,leg2);
	xlabel('Log of Initial Shock, ln \theta_{t}');
	ylabel('\beta E(c_{t+1}^{-\gamma}(\alpha \theta_{t+1}k_{t+1}^{\alpha-1}+1 - \delta)|\theta_t)');
print -noui -deps2 -r300 accuracyt.eps;
end;

diary off;

% **********************************************************************

%       Characteristics of the algorithm:
%
%       1. State variables are the lns of capital,ln(k), and
%          productivity,ln(t).
%       2. The conditional expectation is approximated using
%          exp{ Pn[ln(k),ln(t)]}, where Pn[x] is a Chebyshev polynomial
%          of order n in x.
%       3. Since Chebyshev polynomials are defined on the interval [-1,1]
%          the state variables are scaled to lie in this interval.
%       4. The conditional expectation is calculated numerically using
%          Hermite Gaussian quadrature.
%       5. A functional iteration scheme is used to find the parameters
%          of the approximating function.
%
%       Steps of the algorithm:
%
%       I.   Choose initial values psi0 for the coefficients of the
%            approximating function. Note that for given values of psi
%            it is easy to calculate the values for the consumption and
%            capital choice as a function of the beginning-of-period
%            capital stock and the productivity shock by replacing the
%            true conditional expectation by the approximating function
%            in the Euler equation.
%
%       II.  At each grid point numerically calculate the conditional
%            expectation of
%
%            dfactor*{confun[ln(k'),rho*ln(t)+sigma*epsi']}
%            *{alpha*(k')^(alpha-1)*exp(rho*ln(t)+sigma*innovation')+1-delta}
%
%            where a ' indicates next period's value.
%            Note that this is a function of the random variable epsilon',
%            which is assumed to be a N(0,1) random variable so we can use
%            Hermite Gaussian quadrature. The conditional expectation is
%            also a function of ln(t) which is known at each grid point
%            and ln(k') which is a function of ln(t) and ln(k) and has
%            to be calculated before the conditional expectation is
%            calculated. psi0 is used to calculate next period's capital
%            stock, savfun[ln(k),ln(t)], next period's consumption level
%            for differnt values of the innovation epsi,
%            confun[ln(k'),ln(t')].
%
%       III. Project the ln of the calculated conditional expectations
%            on Pn[ln(k),ln(t)]. Note that this is just a simple
%            linear regression and that taking lns is allowed
%            since there is no stochastic error term in the equation
%            that equates the conditional expectation to the approximating
%            function.
%
%       IV.  A nonlinear equation solver is used to find the values of the
%            coefficients in the approximating function for which the
%            numerically calculated values are equal to approximation. Note
%            that for this program number of knod = po_k + 1 and tnod = po_t+1.
%            Thus, you can get an exact fit at the grid points. You could have
%            more nodes than coefficients but then you have to use a
%            minimization routine.
%
%
%       Functions used:
%
%       1.   confun(lnk1,lnt1,psi): calculates the ln of the consumption level
%            as a function of the unscaled levels of the state variables.
%       2.   savfun(lnk1,lnt1,psi): calculates the ln of the capital stock
%            as a function of the unscaled levels of the state variables.
%       3.   aexp(epsi,lnk2,lnt1,psi):
%            calculates the discounted marginal utility
%            of next period's marginal product of capital as a function of
%            the ln of next period's capital stock (which is known today),
%            the ln of today's productivity level, and next period's innovation
%       4.   scalup(x,xmin,xmax): a linear function to map a variable defined
%            on [-1,1] to [xmin,xmax].
%       5.   scaldown(x,xmin,xmax): a linear function to map a variable defined
%            on [xmin,xmax] to [-1,1].
%       6.   chebnode(n): function to define n Chebyshev nodes.
%       7.   chebpol(ord,p): function to define a (p+1) vector with
%            a constant and the p Chebyshev polynomial terms.
%       8.   makepoly(ord,x): For the (Tx2) matrix x and the (1x2) vector ord
%            this creates a [T x ord(1)*ord(2)] matrix where each row contains
%            the ord(1)*ord(2) bivariate Chebyshev polynomial terms.
%       9.   hernodes(n): function to calculate the (nx2) matrix with the
%            n Hermite nodes and weights.
%       10.  numi(fun,lnk2,lnt1,psi): function to numerically integrate the
%            function fun(epsi) using Hermite Gaussian quadrature when
%            epsi is a normally distributed random variable with zero mean
%            and standard deviation sigma.
%
%============================================================================
%
%      EXAMPLE; Solution for psi when
%
%      po_k = 3, po_t = 3, knod = 7, tnod = 7
%
%      alpha   = 0.33D0
%      dfactor = 0.95D0
%      nu              = 2.0D0
%      delta   = .95D0
%
%      sigma   = .01D0
%      rho             = 0.95D0
%
%
%      kmin    = 0.8D0
%      kmax    = 1.2D0
%
%      tmin    = -3.5D0
%      tmax    =  3.5D0
%
%
%       1.856033636012942
% -2.357064133658491E-001
% -2.439818435654846E-005
%  1.206701028617710E-006
% -1.136367830054171E-001
%  4.039760903712409E-005
% -8.929668254395991E-006
%  1.590293624380350E-007
%  3.307958486399994E-005
%  1.066761973693250E-005
% -2.725911314969676E-007
%  3.537088810880155E-009
% -1.998446325740655E-006
%  2.363115533506463E-007
% -5.022199010994641E-009
%  4.228266668466543E-011
%
%
%

