/* ** peaproi1.prg, GAUSS program to solve the neoclassical growth model ** ** This program needs the procedures in pealib.src to run. ** ** - Projection method & Iterative solution procedure ** ** 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 = 7; /* Gridpoints for k */ tnod = 7; /* 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; /* ---------------------------------------------------------------------- */ /* 3 The Main Loop ** =============== */ psiter = 1; do while psiter <= maxiter; /* 3.1 Compute the expected value (Step II) ** ---------------------------------------- ** Here we calculate the conditional expectation at each grid point. ** First, we calculate next period's capital stock as a function of ** the ln of the current capital stock and the ln of the technology ** shock. Second, we numerically integrate the function aexp over the ** random innovation to the technology shock. psi is an input of the ** numerical integration since it is used to evaluate the consumption ** function at lnk2 and lnt2. */ yvec = zeros(dimgrid,1); i = 1; do while i <= dimgrid; lnk1 = scalup(xgrid[i,1],kmin,kmax); lnt1 = scalup(xgrid[i,2],tmin,tmax); lnk2 = savfun(lnk1,lnt1,psi); yvec[i] = numi(&aexp,lnk2,lnt1,psi); i = i + 1; endo; /* 3.2 Update the Polynomial Coefficients (Step III & IV) ** ------------------------------------------------------ ** Here we calculate new coefficients for the parameterized expectation ** by regressing the ln of yvec(i) on the polynomial terms. Then we ** measure the distance between the new and the old psi and if the ** difference is not small enough then we take a weighted average ** and repeat the loop. */ psin = xr*ln(yvec); delpsi = sqrt((psi-psin)'(psi-psin)); print psiter~delpsi; if delpsi < psitol; break; else; psi = lrate*psin+(1-lrate)*psi; endif; psiter = psiter + 1; /* The following end statement ends the psiter do-loop. */ endo; /* ** 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. If the values for psi that come out of step 3 are close to ** the ones used in step 2 then the algorithm has converged, ** if not step 2 and 3 have to be repeated. ** ** ** 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 ** ** */