/* ** peapro1.prg, GAUSS program to solve the neoclassical growth model ** ** This program needs the procedures in pealib.src to run. ** ** - 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 ** ************************************************************************ */ new; library pgraph; output file=gauss.out reset; /* 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 = .01; /* Standard Deviation for ln noise in technology shock*/ rho = 0.95; /* Persistence of ln 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 ln theta ** in standard deviations away from mean */ tmax = 3.5; /** Highest gridpoint for ln theta ** in standard deviations away from mean */ maxiter = 2000; /* Maximum Iterations to find ** parameters of polynomial */ psitol = 1e-5; /* Convergence criterion */ lrate = 0.4; /* parameter to control updating of ** new projection coefficient ** lrate = 1 means complete updating ** lrate < 1 means partial updating */ /* 1.3 Compute the Steady State ** ---------------------------- */ ks = ( (1-dfactor+delta*dfactor) / (alpha*dfactor) )^(1/(alpha-1)); kmin = ln(ks * kmin); kmax = ln(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 * *************************************************************************** */ #include pealib.src; /* 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); xgrid = zeros(knod*tnod,2); ki = 1; do while ki <= knod; ti = 1; do while ti <= tnod; xgrid[(ki-1)*tnod+ti,.] = kg[ki]~tg[ti]; ti = ti + 1; endo; ki = ki + 1; endo; dimgrid = rows(xgrid); nvars = cols(xgrid); /* 2.2 Create the projection matrix ** -------------------------------- ** In step III of the algorithm we regress the calculated conditional ** expectation terms on the polynomial terms. Since the regressors never ** change we can compute and store the X matrix of the regression and ** the inv(X'X)*X' matrix. */ xmat = makepoly((po_k~po_t),xgrid); xr = inv(xmat'*xmat)*xmat'; /* 2.3 Initialize the coefficients of the polynomial (Step I) ** ----------------------------------------------------------- */ /* 2.3.1. Load the file ** -------------------- */ psisize = (po_k+1)*(po_t+1); psi = zeros(psisize,1); load psi0[psisize,1] = psi0.in; /* 2.3.2. Initialize psi ** --------------------- */ psi = psi0; /* ---------------------------------------------------------------------- */ /* Call the Equation Solver */ {psi,F} = nsecant(&gridfun,psi0,0,xgrid); /* ** 3.3 Print Parameters ** -------------------- */ print "alpha ";; alpha; print "dfactor ";; dfactor; print "nu ";; nu; print "delta ";; delta; print "rho ";; rho; print "sigma ";; sigma; output off; screen off; output file=psi0.out reset; psi0; output off; screen on; output file = gauss.out on; /* ** *************************************************************************** ** * 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 = rndn(simtime+1,1) * sigma; /* 4.1 Start at the Steady State ** ----------------------- */ lnks[1] = ln( ( (1-dfactor+delta*dfactor) / (alpha*dfactor) )^(1/(alpha-1)) ); lnts[1] = 0; print; print "Steady State "; print "Capital: ";;ks[1]; /* 4.1 Run the simulation ** ---------------------- */ ti = 1; do while ti <= 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]; ti = ti + 1; endo; /* 4.2 Evaluate the simulation ** --------------------------- */ print; print "Over the simulation the means and variances were"; print "\t \t Capital \t Consumption \t Shock"; print "Means \t";;meanc(lnks[1:simtime]~lncs~lnts[1:simtime])'; print "Std.D \t";;stdc(lnks[1:simtime]~lncs~lnts[1:simtime])'; print "Min \t";;minc(lnks[1:simtime]~lncs~lnts[1:simtime])'; print "Max \t";;maxc(lnks[1:simtime]~lncs~lnts[1:simtime])'; if graphs == 1; /* figure(1) */ title("Log of Consumption and Capital Stock Over the Simulation"); xlabel("Time"); ylabel("Consumption, Capital"); _plegstr = "Consumption\000Capital Stock"; _plegctl = 1; graphprt("-C=1 -W=20"); xy(1,lncs[1:simtime]~lnks[1:simtime]); dos copy pqgcvt.out simu.eps; graphset; endif; endif; /* This ends the simulation if */ /* 4.3 Check the Accuracy ** ---------------------- */ if accur == 1; tobs = ceil(tnod/2); obs = 50; lnk1 = seqa(kmin,(kmax-kmin)/(obs-1),obs); lnt1 = scalup(tg[tobs]*(ones(obs,1)),tmin,tmax); lnk2 = savfun(lnk1,lnt1,psi); actk = zeros(obs,1); ai = 1; do while ai <= obs; actk[ai,1] = numi(&aexp,lnk2[ai],lnt1[ai],psi); ai = ai + 1; endo; k1ls = scaldown(lnk1,kmin,kmax); t1ls = scaldown(lnt1,tmin,tmax); apoly = makepoly((po_k~po_t),(k1ls~t1ls)); appk = exp(apoly*psi); /* figure(2) */ title("Accuracy Test, Varying Initial Capital Stock"); xlabel("ln of Initial Capital Stock, ln k]t["); _plegstr ="Actual Expectation\000 Approximate Expectation"; _plegctl = 1; graphprt("-C=1 -W=20"); xy(lnk1~lnk1,actk~appk); dos copy pqgcvt.out accuk.eps; graphset; kobs = ceil(knod/2); lnt1 = seqa(tmin,(tmax-tmin)/(obs-1),obs); lnk1 = scalup(kg[kobs]*ones(obs,1),kmin,kmax); lnk2 = savfun(lnk1,lnt1,psi); actt = zeros(obs,1); ai = 1; do while ai <= obs; actt[ai,1] = numi(&aexp,lnk2[ai],lnt1[ai],psi); ai = ai + 1; endo; k1ls = scaldown(lnk1,kmin,kmax); t1ls = scaldown(lnt1,tmin,tmax); apoly = makepoly((po_k~po_t),(k1ls~t1ls)); appt = exp(apoly*psi); /* figure(3) */ fonts("simplex simgrma"); title("Accuracy Test, Varying Productivity Shock"); _plegstr ="Actual Expectation\000 Approximate Expectation"; _plegctl = 1; xlabel("ln of Productivity Shock, ln \202Q\201]t["); graphprt("-C=1 -W=20"); dos copy pqgcvt.out accut.eps; xy(lnt1~lnt1,actt~appt); graphset; endif; output 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 minimizati on ** 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 ** ** */