/*
   file: Uzawa_Lucas.g
   GAUSS code for solving the continuous time Uzawa-Lucas Model without externalities. 
   
   Production functions:
     Goods sector: y = c + kdot + delta*k = A*(k^alpha)*(u*h)^1-alpha
     Education sector: hdot + delta*h = B*(1-u)*h
   Utility function:
     c^(1-theta)-1/(1-theta)  
     
   State-like variable, w = k/h
   Control-like variable, x = c/k
   Control-like variable, u

   The users are advised to refer to Appendix 5B and pp.593-596, Barro, R. and Sala-i-Martin, X. (2004), Economic Growth, MIT Press.
   The program is used at your own risk. If you find it useful in your work, please cite this program.  
   
   Programmed by 
   Cheuk-Yin Ho 
   Department of Economics and Finance 
   City University of Hong Kong
   Email: cheukho@cityu.edu.hk
   Last update: September 2007
*/


clear _par;                               /* for temporary storage */
format /mb1 /rd 8,6;                      /* reset printing format */
output file = Uzawa_Lucas.out on;         /* open output file */

alpha = 0.3;                              /* physical capital income share */
delta = 0.05;                             /* depreciation rate of capital */
rho = 0.04;                               /* consumer's subjective discount rate */
theta = 2;                                /* 1/theta = constant intertemporal elasticity of substitution */ 
A = 1;                                    /* technology level in the goods sector*/
B = 0.13;                                 /* technology level in the education sector*/
par = (alpha|delta|rho|theta|A|B);        /* stack into a  6 x 1 vector  */
q0 = (1 | 1| 1 );                         /* initial guess for the nonlinear eq. solver eqSolve */
qstar = steady(&eqsys,par,q0);            
wstar = qstar[1];                         /* steady-state of w */
xstar = qstar[2];                         /* steady-state of x */
ustar = qstar[3];                         /* steady-state of u */
aa = gradp(&eqsys,qstar);                 /* the matrix of linearized system */
{va,ve} = eigv(aa);                       /* va = eigenvalues; ve = eigenvectors */


{wv, xv, dxv} = exact(qstar, par);        /* wv is the solution of w; xv[.,1] is the solution of x; xv[.,2] is the solution of u */ 
z = A*(xv[1,2]^(1-alpha))*wv[1]^(-(1-alpha));
m = wv[1]*(z-xv[1,1]-B*(1-xv[1,2]));
{t, wt, dg} = rkode(&dwdt, 0, wv[1]|m, 0.1, 700);  /* t = time; wt = time path of w */
timepath = path(wt, wv, xv[.,1], xv[.,2]);
xt = timepath[.,1];                       /* xt = time path of x */ 
ut = timepath[.,2];                       /* ut = time path of u */


/*
   Time paths of the growth rates of different variables 
   w: state-like variable
   x: control-like variable
   u: time allocated to production in the goods sector
   c: consumption
   k: physical capital
   h: human capital
   y: output
   q: broad output
*/

{wgrowth,xgrowth,ugrowth,cgrowth,kgrowth,hgrowth,ygrowth,qgrowth} = tgrowth(par,wt,xt,ut,wstar,xstar,ustar);



@************************ results **************************@

"eigenvalues";
print;
va';
print;
"eigenvectors";
ve;
print; 
"parameters: alpha, delta, rho, theta, A, and B";
print;
_par';
print;
"wstar, xstar, ustar";
print;
wstar~xstar~ustar;
print;
"wt, xt, ut";
print;
wt~xt~ut;
print;
"wgrowth, xgrowth, ugrowth, cgrowth, kgrowth, hgrowth, ygrowth, qgrowth"; 
print;
wgrowth~xgrowth~ugrowth~cgrowth~kgrowth~hgrowth~ygrowth~qgrowth;



@***********************  procedures  *************************@

/* 
   The system of differential equations in Appendix 5B, Barro and Sala-i-Martin (2004). 
   weq = equation (5.64), xeq = equation (5.70), ueq = equation (5.71)
   Inputs: 
              q   3 x 1,  q[1] = w,  q[2] = x,  q[3] = u.
   Outputs: 
              f   3 x 1,  f[1] = weq,  f[2] = xeq,  f[3] = ueq.
*/

proc (1) = eqsys(q);
local alpha, delta, rho, theta, A, B, weq, xeq, ueq, w, x, u, z, f;
alpha = _par[1];                       /* access the parameter values from the global var _par */
delta = _par[2];
rho = _par[3];
theta = _par[4];
A = _par[5];
B = _par[6];
w = q[1];
x = q[2];
u = q[3];
z = A*(u^(1-alpha))*w^(-(1-alpha));
weq = z-x-B*(1-u);
xeq = ((alpha-theta)*z/theta)+x-(1/theta)*(delta*(1-theta)+rho);
ueq = B*((1-alpha)/alpha)+(B*u)-x;
f = ( weq | xeq | ueq );
retp(f);
endp;

/* 
    Find the steady state numerically. 
    Inputs: 
       &f   proc, return the RHS of the system. 
       par  k x 1, parameters for the equation system.
       q0   3 x 1, (w0 | x0 | u0), initial guess for the nonlinear eq solver eqSolve. 
    Output: 
       qstar  3 x 1, (wstar | xstar | ustar), steady state.  
*/ 

proc (1) = steady(&f,par,q0);
local q,retcode;
local f:proc;
_par = par;                          /* store the parameters in _par for proc f */
eqSolveset;
{q, retcode} = eqSolve(&f, q0);
print;
retp( q );
endp;

/*
   4th order Runge-Kutta ODE integrator for k x 1 system dy/dx = f(x,y), x0, y0 given.
   Inputs: 
       &f    proc for evaluating dy/dx = f(x,y), x scalar, y 1 x k, f(x,y) 1 x k.
       x0    scalar, initial x.
       z0    2 x k,  z0[1,.] = y0, z0[2,.] = f(x0,y0).  
       h     scalar, step size. 
       n     scalar, no of solution points. 
   Outputs: 
       x     n x 1, x[1] = x0, x[i] = x[i-1]+h, i > 1. 
       y     n x k, y = y(x) solution.   
       fv    n x k, f(x,y).   
*/

proc (3) = rkode(&f,x0,z0,h,n);
local k,x,y,fv,i,z1,z2,z3,z4;
local f:proc;
k = cols(z0);
x = seqa(x0,h,n);
y = zeros(n,k);
fv = zeros(n,k);
y[1,.] = z0[1,.];
fv[1,.] = z0[2,.];
for i(2,n,1);
  z1 = fv[i-1,.];
  z2 = f(x[i-1]+0.5*h,y[i-1,.]+0.5*h*z1);
  z3 = f(x[i-1]+0.5*h,y[i-1,.]+0.5*h*z2);
  z4 = f(x[i-1]+h,y[i-1,.]+h*z3);
  y[i,.] = y[i-1,.] + h*(z1/6 + z2/3 + z3/3 + z4/6);       /* update y */
  fv[i,.] = f(x[i],y[i,.]);
endfor;
retp( x,y,fv );
endp;

/* 
  We make use of the procedure rkode to find the soultions of w, x, ,and u.
  Outputs: wv is the solution of w; xv[.,1] is the solution of x; xv[.,2] is the solution of u
*/ 

proc (3) = exact(qstar, par); 
local aa, va, ve, a, v, a0,b0, c0, w0, x0, u0, lo, hi, n, h, wv1, wv2, 
xv1, xv2, dxv1, dxv2, dxv, lxv, wv, xv;
aa = gradp(&eqsys,qstar);              /* the matrix of linearized system */
{va,ve} = eigv(aa);                    /* va = k x 1 eigenvalues, ve = k x k eigenvectors */
if va[1] < 0;
   a = va[1];
   v = real(ve[2,1]/ve[1,1])~real(ve[3,1]/ve[1,1]);          /* force it to real */
elseif va[2] < 0;
   a = va[2];
   v = real(ve[2,2]/ve[1,2])~real(ve[3,2]/ve[1,2]);
elseif va[3] < 0;
   a = va[3];
   v = real(ve[2,3]/ve[1,3])~real(ve[3,3]/ve[1,3]);
endif;

w0 = qstar[1];
x0 = qstar[2];
u0 = qstar[3];
lo = w0;
hi = w0*1.9;
n = 500;
h = (hi-lo)/n;
a0 = x0~u0;
{wv1,xv1,dxv1} = rkode(&policy,w0,a0|v,h,n);
hi = w0;
lo = w0*0.1;
h = (hi-lo)/n;
{wv2,xv2,dxv2} = rkode(&policy,w0,a0|v,-h,n);
wv = rev(trimr(wv2,1,0))|wv1;
xv = rev(trimr(xv2,1,0))|xv1;
dxv = rev(trimr(dxv2,1,0))|dxv1;
retp (wv, xv, dxv);
endp;

/*
   Policy functions; f1 = dx/dw and f2 = du/dw
*/

proc (1) = policy(w,y);
local alpha, delta, rho, theta, A, B, weq, xeq, ueq, f1, f2, f, x, u, z;
alpha = _par[1];                      
delta = _par[2];
rho = _par[3];
theta = _par[4];
A = _par[5];
B = _par[6];
x = y[.,1];
u = y[.,2];
z = A*(u^(1-alpha))*w^(-(1-alpha));
weq = z-x-B*(1-u);
xeq = ((alpha-theta)*z/theta)+x-(1/theta)*(delta*(1-theta)+rho);
ueq = B*((1-alpha)/alpha)+(B*u)-x;
f1 = (x*xeq)/(w*weq);
f2 = (u*ueq)/(w*weq);
f = f1~f2;
retp(f);
endp;


/* Policy function; f = dw/dt */

proc(1) = dwdt(t, wr);
local g, j, n, z, alpha, delta, rho, theta, A, B, f, cr, ur, xr;
alpha = _par[1];                      /* access the parameter values from the global var _par */
delta = _par[2];
rho = _par[3];
theta = _par[4];
A = _par[5];
B = _par[6];
g = cols(wr);
xr = zeros(1,g);
ur = zeros(1,g);
j = 1;
do until j > g;
   xr[j] = polyint(wv, xv[.,1], wr[j]);       /* interpolation */
   ur[j] = polyint(wv, xv[.,2], wr[j]);
   j = j+1;
endo;
z = A*(ur^(1-alpha))*wr^(-(1-alpha));
f = wr*(z-xr-B*(1-ur));
retp(f);
endp;


/* Find the time paths of x and u*/

proc(1) = path(wt, wv, xv, uv);
local i,n,q; 
n = rows(wt);
q = zeros(n,2);
i = 1;
do until i > n;
   q[i,1] = polyint(wv, xv, wt[i]);           /* interpolation */
   q[i,2] = polyint(wv, uv, wt[i]);
   i = i+1;
endo;
retp(q);
endp;


/* Find the time paths of wgrowth, xgrowth, ugrowth, cgrowth, kgrowth, hgrowth, ygrowth,and qgrowth */

proc(8) = tgrowth (par,wt,xt,ut, wstar,xstar,ustar);
local n, z, zstar, wgrowth, xgrowth, ugrowth, cgrowth, kgrowth, hgrowth, ygrowth, qgrowth;
alpha = _par[1];                      
delta = _par[2];
rho = _par[3];
theta = _par[4];
A = _par[5];
B = _par[6];
n = rows(wt);
z = zeros(n,1);
zstar = B/alpha;
z[.,1] = A*ut^(1-alpha).*(wt^-(1-alpha));
wgrowth = (z-zstar) - (xt-xstar)+B*(ut-ustar);
xgrowth = ((alpha-theta)/theta)*(z-zstar)+(xt-xstar);
ugrowth = B*(ut-ustar) - (xt-xstar); 
cgrowth = 1/theta*(alpha*z-delta-rho);
kgrowth = cgrowth - xgrowth;
hgrowth = kgrowth - wgrowth;
ygrowth = alpha*kgrowth + (1-alpha)*(ugrowth+hgrowth);
qgrowth = ygrowth - ugrowth.*((1-alpha)/(1-alpha+alpha*ut));
retp(wgrowth,xgrowth,ugrowth,cgrowth,kgrowth,hgrowth,ygrowth,qgrowth);
endp;

output off;
end;