CFPP$ UNROLL (32) F c c SUBROUTINE: CUMOM.F c c TITLE: moment calculation for simulated method of moments c estimation of parameters in cu3da.f model c c PROBLEM: calculate equilibria as in cu3da.f, simulate c time series, calculate certain moments that will c be compared with those in the data c c ALGORITHM: finite element method used to compute equilibria; c simulated method of moments used to estimate parameters c c MODEL: stochastic growth model with government capital c and varying hours and employment; the planner c chooses c1[t],c0[t],i[t],n[t],h[t] to maximize c c E[ sum beta^t {n[t]*u(c1[t],1-h[t])+(1-n[t]) c t *u(c0[t],1)}*(1-a[t]) | x[0]] c c subject to c c c[t]+i[t] = z[t]*f(k[t],kg[t],n[t],h[t])-ig[t]-cg[t] c c[t] = n[t]*c1[t]+(1-n[t])*c0[t] c k[t+1] = (1-delta)*k[t]+i[t] c kg[t+1] = (1-delta)*kg[t]+ig[t] c z[t+1] = z[t]^psi*eps[t] c x[t] = (k[t],kg[t],z[t]) c c FUNCTIONS: u(c,l) = [(c^gamma*l^(1-gamma))^(1-omega)-1]/(1-omega) c F(k,kg,n,h) = lam*z*(phi*k^rho+(1-phi)*kg^rho)^ ... c (theta/rho)*n^(1-theta)*h c c PARAMETERS: param = [beta,delta,gamma,lam,omega,pa,pb,phi,psi, ... c rho,sig,theta]' c c INPUTS: param: initial parameter vector, param c c OUTPUTS: mm: vector of moments c c OTHER (LINPACK) blas.f dgefa.f, dgesl.f, dgeco.f c SUBROUTINES: (KMN) dgefs.f, machcon.f, xerror.f c (McGrattan) qgausl.f c c REFERENCES: Eischen, Jeff, (notes on the Finite Element Method) c Kahaner, Moler, and Nash (1989), NUMERICAL METHODS c AND SOFTWARE (KMN) c Saad, Youcef (1993), ``Lecture notes: Scientific c Computation 8002,'' Univ. Minnesota c c SEE ALSO: cu3da.f,cu3db.f c c Ellen McGrattan, 9-29-93 c Revised, 9-29-93, ERM c subroutine cumom(param,mm) parameter (na=96,nx=15,ny=1,nz=2,ns=5,nas=480,nxm=3,nym=3, c nzm=3,ne1=10,nt=10000,nw=10000) implicit double precision (a-h,o-z) double precision f(nas), df(nas,nas), work(nw), c xa(nx+1), ya(ny+1), za(nz+1), xx(nxm,nx), wx(nxm,nx), c yy(nym,ny), wy(nym,ny), zz(nzm,nz), wz(nzm,nz), c x(nxm), u(nxm), y(nym), v(nym), z(nzm), w(nzm), c c1(nt,ns), b(nt,ns), res(nt,ns), n0(nt,ns), h(nt,ns), c a0(na,ns), a1(na,ns), zzt(ne1), e1(ne1), w1(ne1), c wgt(ne1), sum(9), wrk(ne1*ns,8), mcp(ns,ns), ginv(ns), c cg(ns), army(ns), lam, n0t integer iwork(nw), nxp(nx), nyp(ny), nzp(nz), e(nx*ny*nz,8), c kk(ne1), ind(ne1*ns) c beta = param(1) delta = param(2) gamma = param(3) lam = param(4) omega = param(5) pa = param(6) pb = param(7) phi = param(8) psi = param(9) rho = param(10) sig = param(11) theta = param(12) zeta = param(13) ginv(1) = 0.3e-04 ginv(2) = 0.433e-01 ginv(3) = 0.59e-03 ginv(4) = 0.367e-02 cg(1) = 0.08 cg(2) = 0.19 cg(3) = 0.10 cg(4) = 0.13 army(1) = 0.002 army(2) = 0.0674 army(3) = 0.011 army(4) = 0.020 mcp(1,1) = 0.962 mcp(1,2) = 0.000 mcp(1,3) = 0.019 mcp(1,4) = 0.019 mcp(2,1) = 0.000 mcp(2,2) = 0.660 mcp(2,3) = 0.000 mcp(2,4) = 0.340 mcp(3,1) = 0.043 mcp(3,2) = 0.000 mcp(3,3) = 0.895 mcp(3,4) = 0.062 mcp(4,1) = 0.000 mcp(4,2) = 0.115 mcp(4,3) = 0.400 mcp(4,4) = 0.485 xa(1) = 0.150000E-01 xa(2) = 0.500000E-01 xa(3) = 0.105000E+00 xa(4) = 0.180000E+00 xa(5) = 0.275000E+00 xa(6) = 0.390000E+00 xa(7) = 0.525000E+00 xa(8) = 0.680000E+00 xa(9) = 0.855000E+00 xa(10) = 0.105000E+01 xa(11) = 0.126500E+01 xa(12) = 0.150000E+01 xa(13) = 0.175500E+01 xa(14) = 0.203000E+01 xa(15) = 0.232500E+01 xa(16) = 0.264000E+01 ya(1) = 0.000000E+00 ya(2) = 0.100000E+00 ya(3) = 0.200000E+00 ya(4) = 0.300000E+00 ya(5) = 0.400000E+00 za(1) = 0.900000E+00 za(2) = 0.100000E+01 za(3) = 0.110000E+01 do 10 i=1,nx nxp(i) = 2 10 continue do 20 i=1,ny nyp(i) = 2 20 continue do 30 i=1,nz nzp(i) = 2 30 continue delt1 = 1.d0-delta gamm1 = 1.d0-gamma ome1 = 1.d0-omega phi1 = 1.d0-phi rho1 = 1.d0-rho tdrho = theta/rho thet1 = 1.d0-theta cwgt = gamma*ome1 cwgt1 = cwgt-1.d0 hwgt = gamm1*ome1 hwgt1 = hwgt-1.d0 tiny = 0.000000001d0 pi = 3.14159265d0 crit = 0.000001d0 maxit = 20 nx1 = nx+1 ny1 = ny+1 nz1 = nz+1 c do xx l=1,ns ssg = ig(l)/delta ssk = 0.8 ssc = 0.2 sse = 0.4 ssh = 0.6 do xx t=1,5 d = (sse-pa)*pb dpn = 3.d0*d*(d+dsign(d,0.d0))*pb; zkk = lam*(phi*ssk**rho+phi1*ssg**rho)**tdrho c0 = ssc*(1.d0-ssh)**(hwgt/cwgt1) foc(1) = gamm1/gamma*sse**theta*ssc-zkk*(1.d0-ssh) foc(2) = (1.d0-ssh)**hwgt/ome1-(c0/ssc)**cwgt/ome1-dpn/ c (ssc**cwgt)-gamma*(1.d0-ssh)**hwgt/ssc*(ssc-c0- c zkk*thet1*sse**(-theta)*ssh) foc(3) = sse*ssc+(1.d0-sse)*c0+delta*ssk+ig(l)+cg(l)- c zkk*ssh*sse**thet1 foc(4) = beta*(delt1+theta*phi*ssk**(rho-1.d0)*zkk/ c (phi*ssk**rho+phi1*ssg**rho)*sse**thet1*ssh)-1.d0 do xx i=1,4 dssk = ssk + del(1,i) dssc = ssc + del(2,i) dsse = sse + del(3,i) dssh = ssh + del(4,i) zkk = lam*(phi*dssk**rho+phi1*dssg**rho)**tdrho c0 = dssc*(1.d0-dssh)**(hwgt/cwgt1) d = (dsse-pa)*pb dpn = 3.d0*d*(d+dsign(d,0.d0))*pb foc1p = gamm1/gamma*dsse**theta*dssc-zkk*(1.d0-dssh) foc2p = (1.d0-dssh)**hwgt/ome1-(c0/dssc)**cwgt/ome1-dpn/ (dssc**cwgt)-gamma*(1.d0-dssh)**hwgt/dssc*(dssc- c0-zkk*thet1*dsse**(-theta)*dssh) foc3p = dsse*dssc+(1.d0-dsse)*c0+delta*dssk+ig(l)+cg(l)- zkk*dssh*dsse**thet1 foc4p = beta*(delt1+theta*phi*dssk**(rho-1.d0)*zkk/(phi* dssk**rho+phi1*dssg**rho)*dsse**thet1*dssh)-1.d0 dssk = ssk - del(1,i) dssc = ssc - del(2,i) dsse = sse - del(3,i) dssh = ssh - del(4,i) zkk = lam*(phi*dssk**rho+phi1*dssg**rho)**tdrho c0 = dssc*(1.d0-dssh)**(hwgt/cwgt1) d = (dsse-pa)*pb dpn = 3.d0*d*(d+dsign(d,0.d0))*pb foc1m = gamm1/gamma*dsse**theta*dssc-zkk*(1.d0-dssh) foc2m = (1.d0-dssh)**hwgt/ome1-(c0/dssc)**cwgt/ome1-dpn/ (dssc**cwgt)-gamma*(1.d0-dssh)**hwgt/dssc*(dssc- c0-zkk*thet1*dsse**(-theta)*dssh) foc3m = dsse*dssc+(1.d0-dsse)*c0+delta*dssk+ig(l)+cg(l)- zkk*dssh*dsse**thet1 foc4m = beta*(delt1+theta*phi*dssk**(rho-1.d0)*zkk/(phi* dssk**rho+phi1*dssg**rho)*dsse**thet1*dssh)-1.d0 dfoc(1,i) = 0.5*(foc1p-foc1m)/del(i,i) dfoc(2,i) = 0.5*(foc2p-foc2m)/del(i,i) dfoc(3,i) = 0.5*(foc3p-foc3m)/del(i,i) dfoc(4,i) = 0.5*(foc4p-foc4m)/del(i,i) continue call dgeco(dfoc,4,4,iwork,rcond,work) if (rcond.eq.0.0) c call xerror( 'dgeco error -- singular matrix ',55,-4,0) call dgesl(dfoc,4,4,iwork,foc,0) ssk = ssk - foc(1) ssc = ssc - foc(2) sse = sse - foc(3) ssh = ssh - foc(4) continue ssi = delta*ssk ssz = 0.d0 c c0 = (lam*dexp(ssz)*(phi*ssk**rho+phi1*ssg**rho)**tdrho* c sse**thet1*ssh-ssi-ginv(l)-cg(l)-sse*ssc)/(1.d0-sse) uu = (1.d0-army(l))*(sse*(ssc**cwgt*(1.d0-ssh)**hwgt+ c (1.d0-sse)*c0**cwgt) do xx i=1,6 dssk = ssk + del(1,i) dssg = ssg + del(2,i) dssz = ssz + del(3,i) dssc = ssc + del(5,i) dssi = ssi + del(6,i) dsse = sse + del(7,i) dssh = ssh + del(8,i) c0 = (lam*dexp(dssz)*(phi*dssk**rho+phi1*dssg**rho)**tdrho* c dsse**thet1*dssh-dssi-ginv(l)-cg(l)-dsse*dssc)/ c (1.d0-dsse) up = (1.d0-army(l))*(dsse*(dssc**cwgt*(1.d0-dssh)**hwgt+ c (1.d0-dsse)*c0**cwgt) dssk = ssk - del(1,i) dssg = ssg - del(2,i) dssz = ssz - del(3,i) dssc = ssc - del(5,i) dssi = ssi - del(6,i) dsse = sse - del(7,i) dssh = ssh - del(8,i) c0 = (lam*dexp(dssz)*(phi*dssk**rho+phi1*dssg**rho)**tdrho* c dsse**thet1*dssh-dssi-ginv(l)-cg(l)-dsse*dssc)/ c (1.d0-dsse) um = (1.d0-army(l))*(dsse*(dssc**cwgt*(1.d0-dssh)**hwgt+ c (1.d0-dsse)*c0**cwgt) du(i)= 0.5*(up-um)/del(i,i) do xx j=1,i dssk = ssk + del(1,i) + del(1,j) dssg = ssg + del(2,i) + del(2,j) dssz = ssz + del(3,i) + del(3,j) dssc = ssc + del(5,i) + del(5,j) dssi = ssi + del(6,i) + del(6,j) dsse = sse + del(7,i) + del(7,j) dssh = ssh + del(8,i) + del(8,j) c23456789012345678901234567890123456789012345678901234567890123456789012 c0 = (lam*dexp(dssz)*(phi*dssk**rho+phi1*dssg**rho)**tdrho* c dsse**thet1*dssh-dssi-ginv(l)-cg(l)-dsse*dssc)/ c (1.d0-dsse) upp = (1.d0-army(l))*(dsse*(dssc**cwgt*(1.d0-dssh)**hwgt+ c (1.d0-dsse)*c0**cwgt) dssk = ssk + del(1,i) - del(1,j) dssg = ssg + del(2,i) - del(2,j) dssz = ssz + del(3,i) - del(3,j) dssc = ssc + del(5,i) - del(5,j) dssi = ssi + del(6,i) - del(6,j) dsse = sse + del(7,i) - del(7,j) dssh = ssh + del(8,i) - del(8,j) c0 = (lam*dexp(dssz)*(phi*dssk**rho+phi1*dssg**rho)**tdrho* c dsse**thet1*dssh-dssi-ginv(l)-cg(l)-dsse*dssc)/ c (1.d0-dsse) upm = (1.d0-army(l))*(dsse*(dssc**cwgt*(1.d0-dssh)**hwgt+ c (1.d0-dsse)*c0**cwgt) dssk = ssk - del(1,i) + del(1,j) dssg = ssg - del(2,i) + del(2,j) dssz = ssz - del(3,i) + del(3,j) dssc = ssc - del(5,i) + del(5,j) dssi = ssi - del(6,i) + del(6,j) dsse = sse - del(7,i) + del(7,j) dssh = ssh - del(8,i) + del(8,j) c0 = (lam*dexp(dssz)*(phi*dssk**rho+phi1*dssg**rho)**tdrho* c dsse**thet1*dssh-dssi-ginv(l)-cg(l)-dsse*dssc)/ c (1.d0-dsse) ump = (1.d0-army(l))*(dsse*(dssc**cwgt*(1.d0-dssh)**hwgt+ c (1.d0-dsse)*c0**cwgt) dssk = ssk - del(1,i) - del(1,j) dssg = ssg - del(2,i) - del(2,j) dssz = ssz - del(3,i) - del(3,j) dssc = ssc - del(5,i) - del(5,j) dssi = ssi - del(6,i) - del(6,j) dsse = sse - del(7,i) - del(7,j) dssh = ssh - del(8,i) - del(8,j) c0 = (lam*dexp(dssz)*(phi*dssk**rho+phi1*dssg**rho)**tdrho* c dsse**thet1*dssh-dssi-ginv(l)-cg(l)-dsse*dssc)/ c (1.d0-dsse) umm = (1.d0-army(l))*(dsse*(dssc**cwgt*(1.d0-dssh)**hwgt+ c (1.d0-dsse)*c0**cwgt) du2(i,j) = 0.25*(upp-upm-ump+umm)/(del(i,i)*del(j,j)) du2(j,i) = du2(i,j) continue continue qmat = 0 wmat = 0 rmat = 0 qwr = 0 do xx j=1,8 qwr(1,j) = du2(1,j) qwr(2,j) = du2(2,j) qwr(3,j) = du2(3,j) qwr(4,j) = du2(4,j) qwr(5,j) = du2(5,j) qwr(6,j) = du2(6,j) qwr(7,j) = du2(7,j) qwr(8,j) = du2(8,j) continue do xx i=1,8 work(i) = 0.5*(du(i)-du2(i,1)*ssk-du2(i,2)*ssg-du2(i,4)* c ssc-du2(i,5)*ssi-du2(i,6)*sse-du2(i,7)*ssh) qwr(4,i) = qwr(4,i) + work(i) qwr(i,4) = qwr(i,4) + work(i) continue qwr(4,4) = qwr(4,4) + uu+(work(1)-du(1))*ssk+ c (work(2)-du(2))*ssg+(work(4)-du(4))*ssc+ c (work(5)-du(5))*ssi+(work(6)-du(6))*sse+ c (work(7)-du(7))*ssh do xx i=1,4 do xx j=1,4 qmat(i,j) = qwr(i,j) amat(i,j) = 0.d0 continue do xx i=1,4 do xx j=6,8 wmat(i,j) = qwr(i,j) bmat(i,j) = 0.d0 continue do xx i=6,8 do xx j=6,8 rmat = qwr(i,j) continue sb = dsqrt(beta) amat(1,1) = delt1*sb amat(2,2) = delt1*sb amat(2,3) = ginv(l)*sb amat(3,3) = psi*sb amat(4,4) = sb bmat(1,2) = sb c c Convert to the no-discount, no-cross product case: c do xx j=1,4 ri(1,j) = rmat(1,j) ri(2,j) = rmat(2,j) ri(3,j) = rmat(3,j) ri(4,j) = rmat(4,j) continue call dgedi(ri,4,4,iwork,dummy,01) do xx i=1,4 do xx j=1,4 rw(i,j) = ri(i,1)*wmat(j,1)+ri(i,2)*wmat(j,2)+ c ri(i,3)*wmat(j,3)+ri(i,4)*wmat(j,4) continue do xx i=1,4 do xx j=1,4 qmat(i,j) = qmat(i,j) -wmat(i,1)*rw(1,j)-wmat(i,2)* c rw(2,j)-wmat(i,3)*rw(j,3)-wmat(i,4)*rw(j,4) continue c c Doubling algorithm used to solve fixed point of Riccati equation c nsv = 4 ndv = 4 A*B*C lxm mxn nxp do xx i=1,nsv do xx j=1,ndv sum1 = 0.d0 do k=1,ndv sum1 = sum1 + bmat(i,k)*ri(k,j) continue wmat(i,j) = sum1 continue do xx i=1,nsv do xx j=1,nsv sum1 = 0.d0 do k=1,ndv sum1 = sum1 + wmat(i,k)*bmat(j,k) continue b0(i,j) = sum1 continue do xx i=1,ndv do xx j=1,nsv sum1 = 0.d0 do k=1,nsv sum1 = sum1 + bmat(k,i)*qmat(k,j) continue wmat(j,i) = sum1 continue do xx i=1,ndv do xx j=1,ndv sum1 = 0.d0 do k=1,nsv sum1 = sum1 + wmat(k,i)*bmat(k,j) continue rb(i,j) = sum1+rmat(i,j) continue call dgeco(rb,ndv,ndv,iwork,rcond,work) if (rcond.eq.0.0) c call xerror( 'dgeco error -- singular matrix ',55,-4,0) call dgesl(rb,ndv,ndv,f0, ...) b0 = bmat*ri*bmat'; rb = bmat'*qmat*bmat+rmat f0 = bmat'*qmat*amat d1 = 1; it = 1; while (d1>tol & it<=maxit); v1 = eye(m)+b0*qmat call dgedi(v1,4,4,iwork,dummy,01) amat= amat*v1*amat; b0 = b0+amat*v1*b0*amat'; qmat = qmat+amat'*qmat*v1*amat; f1 = bmat'*qmat*amat rb = bmat'*qmat*bmat+rmat call dgesl(rb,4,4,f1, ...) d1 = norm(f1-f0); it = it+1; f0 = f1; end rule = -(f1+rw); t = 1 do xx k=1,nz1 do xx j=1,ny1 do xx i=1,nx1 a0(t,l) = rule(1,1)*xa(i)+rule(1,2)*ya(j)+ c rule(1,3)*za(k)+rule(1,4) t = t+1 continue continue -- last l c sum1 = 0.d0 sum2 = 0.d0 call qgausl(ne1,-2.88d0*sig,2.88d0*sig,e1,w1) do 60 i=1,ne1 wgt(i) = dexp(-0.5d0*e1(i)*e1(i)/(sig*sig))/ c dsqrt(2.d0*pi*sig*sig)*w1(i) sum1 = sum1 + wgt(i) 60 continue do 70 i=1,ne1 wgt(i) = wgt(i)/sum1 sum2 = sum2 + wgt(i) 70 continue n = 0 do 80 k=1,nz do 80 j=1,ny do 80 i=1,nx ll = (k-1)*nx*ny+(j-1)*nx+i e(ll,1) = (k-1)*nx1*ny1+(j-1)*nx1+i e(ll,2) = i e(ll,3) = j e(ll,4) = k e(ll,5) = nxp(i) e(ll,6) = nyp(j) e(ll,7) = nzp(k) e(ll,8) = n n = n+nxp(i)*nyp(j)*nzp(k) 80 continue do 90 i=1,nx call qgausl(nxp(i),xa(i),xa(i+1),x,u) do 90 j=1,nxp(i) xx(j,i) = x(j) wx(j,i) = u(j) 90 continue do 100 i=1,ny call qgausl(nyp(i),ya(i),ya(i+1),y,v) do 100 j=1,nyp(i) yy(j,i) = y(j) wy(j,i) = v(j) 100 continue do 110 i=1,nz call qgausl(nzp(i),za(i),za(i+1),z,w) do 110 j=1,nzp(i) zz(j,i) = z(j) wz(j,i) = w(j) 110 continue ne = nx*ny*nz dd = 1.d0 do 120 it=1,maxit do 130 i=1,nas f(i) = 0.d0 do 130 j=1,nas df(i,j) = 0.d0 130 continue do 140 n=1,ne n1 = e(n,5) n2 = e(n,6) n3 = e(n,7) i1 = e(n,8) x1 = xa(e(n,2)) x2 = xa(e(n,2)+1) dx = x2-x1 y1 = ya(e(n,3)) y2 = ya(e(n,3)+1) dy = y2-y1 z1 = za(e(n,4)) z2 = za(e(n,4)+1) dz = z2-z1 nd1= e(n,1) nd2= nd1+1 nd4= nd1+nx1 nd3= nd4+1 nd5= nd1+nx1*ny1 nd6= nd5+1 nd8= nd5+nx1 nd7= nd8+1 do 150 j=1,n2 y(j) = yy(j,e(n,3)) 150 continue do 160 k=1,n3 z(k) = zz(k,e(n,4)) 160 continue c create vectors a,b,c for consumption, leisure, next capital do 170 i=1,n1 x(i) = xx(i,e(n,2)) xp = (x2 - x(i))/dx xm = (x(i) - x1)/dx do 170 j=1,n2 yp = (y2 - y(j))/dy ym = (y(j) - y1)/dy do 170 k=1,n3 zp = (z2 - z(k))/dz zm = (z(k) - z1)/dz ll = i1+(k-1)*n1*n2+(j-1)*n1+i do 170 l=1,ns c1(ll,l) = a0(nd1,l)*xp*yp*zp+a0(nd2,l)*xm*yp*zp+ c a0(nd3,l)*xm*ym*zp+a0(nd4,l)*xp*ym*zp+ c a0(nd5,l)*xp*yp*zm+a0(nd6,l)*xm*yp*zm+ c a0(nd7,l)*xm*ym*zm+a0(nd8,l)*xp*ym*zm c1(ll,l) = dmax1(c1(ll,l),tiny) b(ll,l) = z(k)*lam*(phi*x(i)**rho+(1.d0-phi)* c y(j)**rho)**(theta/rho) 170 continue nn = n1*n2*n3 do 180 ll=i1+1,i1+nn do 180 l=1,ns hrs = 0.4d0 emp = 0.2d0 do 190 m=1,5 c0 = c1(ll,l)*(1.d0-hrs)**(hwgt/cwgt1) c pn = abs((emp*pb-pa*pb)**3.0)+(emp*pb-pa*pb)**3.0 d = (emp-pa)*pb dpn = 3.d0*d*(d+dsign(d,0.d0))*pb d2pn = 6.d0*(d+dsign(d,0.d0))*pb*pb fj11 = -d2pn/(c1(ll,l)**cwgt)-gamma*(1.d0-hrs)** c hwgt/c1(ll,l)*b(ll,l)*thet1*theta*emp** c (-theta-1.d0)*hrs fj22 = b(ll,l) fj12 = hwgt*(1.d0-hrs)**hwgt1/ome1-gamma*hwgt* c (1.d0-hrs)**hwgt1/c1(ll,l)*(c1(ll,l)-c0- c b(ll,l)*thet1*emp**(-theta)*hrs)-gamma* c (1.d0-hrs)**hwgt/c1(ll,l)*b(ll,l)*thet1* c emp**(-theta) fj21 = -gamm1*theta*emp**(-thet1)*c1(ll,l)/gamma det = fj11*fj22-fj12*fj21 f1 = gamm1/gamma*emp**theta*c1(ll,l)-b(ll,l)* c (1.d0-hrs) f2 = (1.d0-hrs)**hwgt/ome1-(c0/c1(ll,l))**cwgt/ c ome1-dpn/(c1(ll,l)**cwgt)-gamma*(1.d0-hrs)** c hwgt/c1(ll,l)*(c1(ll,l)-c0-b(ll,l)*thet1* c emp**(-theta)*hrs) hrs = hrs - (fj11*f1+fj21*f2)/det emp = emp - (fj12*f1+fj22*f2)/det hrs = dmax1(hrs,tiny) emp = dmax1(emp,tiny) 190 continue h(ll,l) = dmin1(hrs,1.d0) n0(ll,l) = dmin1(emp,1.d0) 180 continue 140 continue do 200 l=1,ns do 210 n=1,ne n1 = e(n,5) n2 = e(n,6) n3 = e(n,7) i1 = e(n,8) x1 = xa(e(n,2)) x2 = xa(e(n,2)+1) dx = x2-x1 y1 = ya(e(n,3)) y2 = ya(e(n,3)+1) dy = y2-y1 z1 = za(e(n,4)) z2 = za(e(n,4)+1) dz = z2-z1 nd1= e(n,1) nd2= nd1+1 nd4= nd1+nx1 nd3= nd4+1 nd5= nd1+nx1*ny1 nd6= nd5+1 nd8= nd5+nx1 nd7= nd8+1 do 220 i=1,n1 x(i) = xx(i,e(n,2)) u(i) = wx(i,e(n,2)) 220 continue do 230 j=1,n2 y(j) = yy(j,e(n,3)) v(j) = wy(j,e(n,3)) 230 continue do 240 k=1,n3 z(k) = zz(k,e(n,4)) w(k) = wz(k,e(n,4)) zp = (z2 - z(k))/dz zm = (z(k) - z1)/dz do 240 j=1,n2 yp = (y2 - y(j))/dy ym = (y(j) - y1)/dy do 240 i=1,n1 ll = (k-1)*n1*n2+(j-1)*n1+i+i1 xp = (x2 - x(i))/dx xm = (x(i) - x1)/dx c0 = c1(ll,l)*(1.d0-h(ll,l))**(hwgt/cwgt1) cg1 = cg(l)+zeta*dlog(z(k)) xt = b(ll,l)*n0(ll,l)**thet1*h(ll,l)+delt1* c x(i)-ginv(l)-cg1-n0(ll,l)*c1(ll,l)- c (1.d0-n0(ll,l))*c0 yt = (1-delta)*y(j)+ginv(l) ii = 1 do 250 m=1,nx j1 = int( dmax1(xt-xa(m),0.d0)/(xt-xa(m))+.2 ) ii = ii + j1 250 continue ii = max(1,ii-1) jj = 1 do 260 m=1,ny j1 = int( dmax1(yt-ya(m),0.d0)/(yt-ya(m))+.2) jj = jj + j1 260 continue jj = max(1,jj-1) CFPP$ UNROLL do 270 m=1,ne1 zzt(m) = dexp(psi*dlog(z(k))+e1(m)) kk(m) = 1 do 280 mm=1,nz j1 = int( dmax1(zzt(m)-za(mm),0.d0)/ c (zzt(m)-za(mm))+.2) kk(m) = kk(m) + j1 280 continue kk(m) = max(1,kk(m)-1) 270 continue do 290 m=1,9 sum(m) = 0.d0 290 continue do 300 m=1,nas work(m) = 0.d0 300 continue xy = phi*xt**rho+(1.d0-phi)*yt**rho xr = xt**(rho-1.d0) h1 = 1.d0-h(ll,l) crat = c0/c1(ll,l) d = (n0(ll,l)-pa)*pb dpn = 3.d0*d*(d+dsign(d,0.d0))*pb d2pn = 6.d0*(d+dsign(d,0.d0))*pb*pb tem1 = gamm1/h1*(crat**cwgt-thet1*h1**hwgt1+ c ome1*dpn/(c1(ll,l)**cwgt)-d2pn*n0(ll,l)/ c (gamm1*theta*c1(ll,l)**cwgt)) tem2 = (d2pn*n0(ll,l)/theta+dpn*cwgt1)/ c (c1(ll,l)**(cwgt+1.d0)) dha = tem2/tem1 tem1 = c1(ll,l)/h1*(-gamm1*thet1*h(ll,l)*n0(ll,l)/ c (gamma*theta*h1)+n0(ll,l)/theta*(1.d0-crat)+ c gamm1*n0(ll,l)/gamma+crat*(1.d0-n0(ll,l))* c hwgt/cwgt1) tem2 = -gamm1*thet1*h(ll,l)*n0(ll,l)/(gamma*theta* c h1)+thet1*n0(ll,l)*(1.d0-crat)/theta-crat dxa = tem1*dha+tem2 do 310 mm=1,ne1 do 310 m=1,ns nct = (m-1)*ne1+mm zt = zzt(mm) nn = (kk(mm)-1)*nx*ny+(jj-1)*nx+ii x1t = xa(e(nn,2)) x2t = xa(e(nn,2)+1) dxt = x2t-x1t y1t = ya(e(nn,3)) y2t = ya(e(nn,3)+1) dyt = y2t-y1t z1t = za(e(nn,4)) z2t = za(e(nn,4)+1) dzt = z2t-z1t xpt = (x2t-xt)/dxt xmt = (xt-x1t)/dxt ypt = (y2t-yt)/dyt ymt = (yt-y1t)/dyt zpt = (z2t-zt)/dzt zmt = (zt-z1t)/dzt nd1t = e(nn,1) nd2t = nd1t+1 nd4t = nd1t+nx1 nd3t = nd4t+1 nd5t = nd1t+nx1*ny1 nd6t = nd5t+1 nd8t = nd5t+nx1 nd7t = nd8t+1 ind(nct) = nd1t c1t = a0(nd1t,m)*xpt*ypt*zpt+ c a0(nd2t,m)*xmt*ypt*zpt+ c a0(nd3t,m)*xmt*ymt*zpt+ c a0(nd4t,m)*xpt*ymt*zpt+ c a0(nd5t,m)*xpt*ypt*zmt+ c a0(nd6t,m)*xmt*ypt*zmt+ c a0(nd7t,m)*xmt*ymt*zmt+ c a0(nd8t,m)*xpt*ymt*zmt c1t = dmax1(c1t,tiny) bt = zt*lam*(phi*xt**rho+(1.d0-phi)* c yt**rho)**(theta/rho) hrs = 0.4d0 emp = 0.2d0 do 320 mt=1,5 c0t = c1t*(1.d0-hrs)**(hwgt/cwgt1) d = (emp-pa)*pb dpn = 3.d0*d*(d+dsign(d,0.d0))*pb d2pn = 6.d0*(d+dsign(d,0.d0))*pb*pb fj11 = -d2pn/(c1t**cwgt)-gamma*(1.d0-hrs)** c hwgt/c1t*bt*thet1*theta*emp**(-theta- c 1.d0)*hrs fj22 = bt dc0 = -c1t*(hwgt/cwgt1)*(1.d0-hrs)**(hwgt/ c cwgt1-1.d0) fj12 = hwgt*(1.d0-hrs)**hwgt1/ome1-gamma*hwgt* c (1.d0-hrs)**hwgt1/c1t*(c1t-c0t-bt* c thet1*emp**(-theta)*hrs)-gamma*(1.d0- c hrs)**hwgt/c1t*bt*thet1*emp**(-theta) fj21 = -gamm1*theta*emp**(-thet1)*c1t/gamma det = fj11*fj22-fj12*fj21 f1 = gamm1/gamma*emp**theta*c1t-bt* c (1.d0-hrs) f2 = (1.d0-hrs)**hwgt/ome1-(c0t/c1t)**cwgt/ c ome1-dpn/(c1t**cwgt)-gamma*(1.d0- c hrs)**hwgt/c1t*(c1t-c0t-bt*thet1*emp** c (-theta)*hrs) hrs = hrs - (fj11*f1+fj21*f2)/det emp = emp - (fj12*f1+fj22*f2)/det hrs = dmax1(hrs,tiny) emp = dmax1(emp,tiny) 320 continue ht = dmin1(hrs,1.d0) n0t = dmin1(emp,1.d0) d = (n0t-pa)*pb dpnt = 3.d0*d*(d+dsign(d,0.d0))*pb d2pnt = 6.d0*(d+dsign(d,0.d0))*pb*pb dc1t = (-a0(nd1t,m)*ypt*zpt+a0(nd2t,m)*ypt*zpt c +a0(nd3t,m)*ymt*zpt-a0(nd4t,m)*ymt*zpt c -a0(nd5t,m)*ypt*zmt+a0(nd6t,m)*ypt*zmt c +a0(nd7t,m)*ymt*zmt-a0(nd8t,m)*ymt*zmt) c /dxt out = beta*zt*lam*xy**(theta/rho)*n0t**thet1*ht trm1 = beta*delt1+out/xy*phi*theta*xr h1 = 1.d0-ht crat = c0t/c1t tem1 = gamm1/h1*(crat**cwgt-thet1*h1**hwgt1+ c ome1*dpnt/(c1t**cwgt)-d2pnt*n0t/ c (gamm1*theta*c1t**cwgt)) tem2 = (d2pnt*n0t/theta+dpnt*cwgt1)/(c1t** c (cwgt+1.d0)) tem3 =-d2pnt*phi*xr*n0t/(xy*c1t**cwgt) dhta = (tem2*dc1t+tem3)*dxa/tem1 dnta =-n0t*dhta/(theta*h1)-n0t*dc1t*dxa/(theta* c c1t)+phi*xr*n0t*dxa/xy dtm1a = phi*theta*xr*out/xy*((-rho1/xt+(theta/ c rho-1.d0)*phi*rho*xr/xy)*dxa+thet1/n0t* c dnta+dhta/ht) dhtb = tem2/tem1 dntb =-n0t*(dhtb/h1+1.d0/c1t)/theta dtm1b = phi*theta*xr*out/xy*(thet1/n0t*dntb+ c dhtb/ht) trm2 = gamma*c1t**cwgt1*h1**hwgt*(1.d0-army(m)) dtm2a = trm2*(cwgt1/c1t*dc1t*dxa-hwgt/h1*dhta) dtm2b = trm2*(cwgt1/c1t-hwgt/h1*dhtb) c sum(1) = sum(1) + mcp(l,m)*wgt(mm)*trm1*trm2 c drhs = mcp(l,m)*wgt(mm)*(dtm1a*trm2+trm1*dtm2a) sum(2) = sum(2) + drhs*xp*yp*zp sum(3) = sum(3) + drhs*xm*yp*zp sum(4) = sum(4) + drhs*xm*ym*zp sum(5) = sum(5) + drhs*xp*ym*zp sum(6) = sum(6) + drhs*xp*yp*zm sum(7) = sum(7) + drhs*xm*yp*zm sum(8) = sum(8) + drhs*xm*ym*zm sum(9) = sum(9) + drhs*xp*ym*zm c drhs = mcp(l,m)*wgt(mm)*(dtm1b*trm2+trm1*dtm2b) wrk(nct,1) = drhs*xpt*ypt*zpt wrk(nct,2) = drhs*xmt*ypt*zpt wrk(nct,3) = drhs*xmt*ymt*zpt wrk(nct,4) = drhs*xpt*ymt*zpt wrk(nct,5) = drhs*xpt*ypt*zmt wrk(nct,6) = drhs*xmt*ypt*zmt wrk(nct,7) = drhs*xmt*ymt*zmt wrk(nct,8) = drhs*xpt*ymt*zmt 310 continue dlhs = -gamma*c1(ll,l)**cwgt1*(1.d0-h(ll,l))**hwgt* c (1.d0-army(l))*(cwgt1/c1(ll,l)-hwgt/ c (1.d0-h(ll,l))*dha) m = (l-1)*na work(m+nd1) = sum(2)+dlhs*xp*yp*zp work(m+nd2) = sum(3)+dlhs*xm*yp*zp work(m+nd3) = sum(4)+dlhs*xm*ym*zp work(m+nd4) = sum(5)+dlhs*xp*ym*zp work(m+nd5) = sum(6)+dlhs*xp*yp*zm work(m+nd6) = sum(7)+dlhs*xm*yp*zm work(m+nd7) = sum(8)+dlhs*xm*ym*zm work(m+nd8) = sum(9)+dlhs*xp*ym*zm do 330 mm=1,ne1 do 330 m=1,ns nct = (m-1)*ne1+mm nd1t = ind(nct) nd2t = nd1t+1 nd4t = nd1t+nx1 nd3t = nd4t+1 nd5t = nd1t+nx1*ny1 nd6t = nd5t+1 nd8t = nd5t+nx1 nd7t = nd8t+1 mt = (m-1)*na work(mt+nd1t) = work(mt+nd1t) + wrk(nct,1) work(mt+nd2t) = work(mt+nd2t) + wrk(nct,2) work(mt+nd3t) = work(mt+nd3t) + wrk(nct,3) work(mt+nd4t) = work(mt+nd4t) + wrk(nct,4) work(mt+nd5t) = work(mt+nd5t) + wrk(nct,5) work(mt+nd6t) = work(mt+nd6t) + wrk(nct,6) work(mt+nd7t) = work(mt+nd7t) + wrk(nct,7) work(mt+nd8t) = work(mt+nd8t) + wrk(nct,8) 330 continue res(ll,l)= (-gamma*c1(ll,l)**cwgt1*(1.d0-h(ll,l))** c hwgt*(1.d0-army(l))+sum(1))*u(i)*v(j)*w(k) m = (l-1)*na f(m+nd1) = f(m+nd1) + res(ll,l)*xp*yp*zp f(m+nd2) = f(m+nd2) + res(ll,l)*xm*yp*zp f(m+nd3) = f(m+nd3) + res(ll,l)*xm*ym*zp f(m+nd4) = f(m+nd4) + res(ll,l)*xp*ym*zp f(m+nd5) = f(m+nd5) + res(ll,l)*xp*yp*zm f(m+nd6) = f(m+nd6) + res(ll,l)*xm*yp*zm f(m+nd7) = f(m+nd7) + res(ll,l)*xm*ym*zm f(m+nd8) = f(m+nd8) + res(ll,l)*xp*ym*zm c tem1 = u(i)*v(j)*w(k) do 340 mm=1,nas tem2 = tem1*work(mm) df(m+nd1,mm) = df(m+nd1,mm) + tem2*xp*yp*zp df(m+nd2,mm) = df(m+nd2,mm) + tem2*xm*yp*zp df(m+nd3,mm) = df(m+nd3,mm) + tem2*xm*ym*zp df(m+nd4,mm) = df(m+nd4,mm) + tem2*xp*ym*zp df(m+nd5,mm) = df(m+nd5,mm) + tem2*xp*yp*zm df(m+nd6,mm) = df(m+nd6,mm) + tem2*xm*yp*zm df(m+nd7,mm) = df(m+nd7,mm) + tem2*xm*ym*zm df(m+nd8,mm) = df(m+nd8,mm) + tem2*xp*ym*zm 340 continue 240 continue 210 continue 200 continue itask = 1 call dgefs (df, nas, nas, f, itask, ind1, work, iwork, rcond) sum3 = 0.d0 do 350 j=1,ns do 350 i=1,na tem1 = f((j-1)*na+i) a1(i,j) = a0(i,j) - tem1 sum3 = sum3 + tem1*tem1 350 continue dd = dsqrt(sum3)/float(nas) if (dd.lt.crit) go to 999 do 360 j=1,ns do 360 i=1,na a0(i,j) = a1(i,j) 360 continue 120 continue 999 continue --> calculate moments % for each state, compute the hours (h), the fraction working (n0), % and consumption for those not working (c0) xx = kron(ones(ny1*nz1,1),x)*ones(1,ns); yy = kron(ones(nz1,1),kron(y,ones(nx1,1)))*ones(1,ns); zz = kron(z,ones(nx1*ny1,1))*ones(1,ns); ii = ones(n,1)*ig'; bb = ones(n,1)*cg'; aa = ones(n,1)*army'; h = 0.5*ones(n,ns); n0 = 0.1*ones(n,ns); zkk = lam*zz.*(phi*xx.^rho+(1-phi)*yy.^rho).^(theta/rho); for m=1:30 c0 = c1.*(1-h).^(hwgt/cwgt1); tem = (n0-pa)*pb; dpn = 3*tem.*tem.*(1+sign(tem))*pb; d2pn = 6*tem.*(1+sign(tem))*pb*pb; fj11 =-d2pn./(c1.^cwgt)-gamma*thet1*theta*(1-h).^hwgt./c1.*zkk.*n0.^ ... (-theta-1).*h; fj22 = zkk; fj12 = hwgt/ome1.*(1-h).^hwgt1-gamma*hwgt*(1-h).^hwgt1./c1.*(c1-c0- ... thet1*zkk.*n0.^(-theta).*h)-gamma*thet1*(1-h).^hwgt./c1.*zkk.* ... n0.^(-theta); fj21 = -gamm1*theta*n0.^(-thet1).*c1/gamma; det = fj11.*fj22-fj12.*fj21; f1 = gamm1/gamma*n0.^theta.*c1-zkk.*(1-h); f2 = (1-h).^hwgt/ome1-(c0./c1).^cwgt/ome1-dpn./(c1.^cwgt)-gamma* ... (1-h).^hwgt./c1.*(c1-c0-thet1*zkk.*n0.^(-theta).*h); h = h - (fj11.*f1+fj21.*f2)./det; n0 = n0 - (fj12.*f1+fj22.*f2)./det; end; % steady states zss = 1; igss = ig(1); cgss = cg(1); ass = army(1); kgss = igss/delta; par = [X(1:12);igss;cgss;ass]; ss = secant('sscu',[.6;.4;.2;.8],par); hss = ss(1); nss = ss(2); c1ss = ss(3); kss = ss(4); % simulations: T = 1000; xt(1) = kss; yt(1) = kgss; zt(1) = zss; igt(1) = igss; cgt(1) = cgss; at(1) = ass; randn('seed',0); et = sig*randn(T,1); st = markov(pi,T); st = st(:); nsim = 1; %load simil %[T,nsim]= size(simil); %st = simil; xt = xt*ones(1,nsim); yt = yt*ones(1,nsim); zt = zt*ones(1,nsim); igt = igt*ones(1,nsim); cgt = cgt*ones(1,nsim); at = at*ones(1,nsim); for j=1:nsim; for t=1:T; l = st(t,j); igt(t,j) = ig(l); cgt(t,j) = cg(l)+zeta*log(zt(t,j)); at(t,j) = army(l); low = find(xt(t,j)>x); high = find(xt(t,j)y); high = find(yt(t,j)z); high = find(zt(t,j)