CFPP$ UNROLL (32) F c c SUBROUTINE: CUMCMOM.F c c TITLE: moment calculation for simulated method of moments c estimation of parameters in cumc3da.f model c c PROBLEM: calculate equilibria as in cumc3da.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 varying hours and employment, and moving costs; c the planner chooses c1[t],c0[t],i[t],n[t],h[t] c 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)-P(n[t])}*(1-a[t])|x[0]] c c subject to c c c[t]+i[t]+m[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 m[t] = M(n[t],n[t-1]) 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],n[t-1],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*(k+kg)^theta*n^(1-theta)*h c M(n,e) = alpha*(n-e)*(n-e) c P(n) = [(n-pa)*pb]^3+abs((n-pa)*pb)^3 c c PARAMETERS: param = [alpha,beta,delta,gamma,lam,omega,pa,pb, ... c psi,sig,theta,zeta]' c c INPUTS: param: initial parameter vector, param c c OUTPUTS: mm: vector of moments c c OTHER (LINPACK) blas.f, dgeco.f, dgedi.f dgefa.f, dgesl.f c SUBROUTINES: (McGrattan) qgausl.f c c REFERENCES: Eischen, Jeff, (notes on the Finite Element Method) c Saad, Youcef (1993), ``Lecture notes: Scientific c Computation 8002,'' Univ. Minnesota c c SEE ALSO: cumc3da.f, cumc3db.f c c Ellen McGrattan, 9-29-93 c Revised, 10-2-93, ERM c subroutine cumcmom(npar,param,obj,grad) parameter (na=63,nx=6,ny=2,nz=2,ns=4,nas=504,nxm=3,nym=3,nzm=3, c ne1=5,nw=20000,nsv=4,ndv=4,nv=8,nt=1140,nm=13,np=11) 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 a0(na,ns), a1(na,ns), b0(na,ns), b1(na,ns), zzt(ne1), c e1(ne1), w1(ne1), wgt(ne1), sum(34), wrk(ne1*ns,32), c mcp(ns,ns), ginv(ns), cg(ns), army(ns), lam, n0, n0t, c mu, wrk2(nas), param(npar), foc(4), dfoc(4,4), del(8,8), c mom, amat(nsv,nsv), bmat(nsv,ndv), qmat(nsv,nsv), c rmat(ndv,ndv), wmat(nsv,ndv), aa(nsv,nsv), bb(nsv,nsv), c cc(nsv,nsv), dd(nsv,nsv), ee(nsv,nsv), qwr(nv,nv), c f0(ndv,nsv), f1(ndv,nsv), du(8), du2(8,8), ri(ndv,ndv), c rw(ndv,nsv), rb(ndv,ndv), det(2), eps(nt), hdat(nm), c hmod(nm), wt(nm,nm), grad(npar), delp(nm,nm+1), st(ns), c outl(57),capl(57),empl(57),hrsl(57),tshl(57), c hmhld(nm), grse(nm,np), ff(nm,np), gg(np,np), se(np) integer iwork(nw), nxp(nx), nyp(ny), nzp(nz), kk(ne1), c ind(ne1*ns), istat(nt) c open(unit=5, file='shocks.dat') open(unit=7, file='tech.dat') do 870 t=1,nt read(5,*) eps(t),istat(t) 870 continue c beta = 0.96 delta = 0.0793 lam = 1.0 omega = 1.001000000000000 isterr = 0 c alpha = param(1) gamma = param(2) pa = param(3) pb = param(4) psi = param(5) sig = param(6) theta = param(7) zeta = param(8) cg(2) = param(9) cg(3) = param(10) cg(4) = param(11) c cg(1) = cg(2)*0.025/0.1872 c cg(5) = cg(4) c army(5) = 0.03187 ginv(1) = 0.0 ginv(2) = 0.0 ginv(3) = 0.0 ginv(4) = 0.0 c ginv(5) = 0.0 c wt(1,1) = 0.0365875 wt(1,2) = -0.0197321 wt(1,3) = 0.00721125 wt(1,4) = -0.00558969 wt(1,5) = -0.0275395 wt(1,6) = 0.0486932 wt(1,7) = -0.0200748 wt(1,8) = 0.0107957 wt(1,9) = -0.00311909 wt(1,10) = -0.00277854 wt(1,11) = 0.00113073 wt(1,12) = -0.000246653 wt(1,13) = 0.000106105 wt(2,1) = -0.0197321 wt(2,2) = 0.0172928 wt(2,3) = -0.00526526 wt(2,4) = 0.00113878 wt(2,5) = 0.024166 wt(2,6) = -0.0271507 wt(2,7) = 0.00255616 wt(2,8) = -0.00703725 wt(2,9) = 0.00241742 wt(2,10) = 0.00216669 wt(2,11) = -0.000439147 wt(2,12) = 0.000225557 wt(2,13) = 0.000102069 wt(3,1) = 0.00721125 wt(3,2) = -0.00526526 wt(3,3) = 0.00513038 wt(3,4) = -0.00385749 wt(3,5) = -0.00737832 wt(3,6) = 0.0212138 wt(3,7) = -0.014193 wt(3,8) = 0.00224916 wt(3,9) = -0.000777001 wt(3,10) = -0.000694728 wt(3,11) = -0.000375799 wt(3,12) = -0.000108527 wt(3,13) = 2.63266e-05 wt(4,1) = -0.00558969 wt(4,2) = 0.00113878 wt(4,3) = -0.00385749 wt(4,4) = 0.00547057 wt(4,5) = 0.00157852 wt(4,6) = -0.0198746 wt(4,7) = 0.0186035 wt(4,8) = -0.00214303 wt(4,9) = 0.000390261 wt(4,10) = 0.000344316 wt(4,11) = 0.000317693 wt(4,12) = 1.43079e-05 wt(4,13) = -7.82474e-05 wt(5,1) = -0.0275395 wt(5,2) = 0.024166 wt(5,3) = -0.00737832 wt(5,4) = 0.00157852 wt(5,5) = 0.033785 wt(5,6) = -0.0379809 wt(5,7) = 0.00358926 wt(5,8) = -0.0098033 wt(5,9) = 0.00337597 wt(5,10) = 0.00302563 wt(5,11) = -0.000621683 wt(5,12) = 0.000316808 wt(5,13) = 0.000142858 wt(6,1) = 0.0486932 wt(6,2) = -0.0271507 wt(6,3) = 0.0212138 wt(6,4) = -0.0198746 wt(6,5) = -0.0379809 wt(6,6) = 0.110596 wt(6,7) = -0.0724617 wt(6,8) = 0.0155384 wt(6,9) = -0.00427207 wt(6,10) = -0.0038096 wt(6,11) = 0.000161267 wt(6,12) = -0.00045724 wt(6,13) = 0.000118812 wt(7,1) = -0.0200748 wt(7,2) = 0.00255616 wt(7,3) = -0.014193 wt(7,4) = 0.0186035 wt(7,5) = 0.00358926 wt(7,6) = -0.0724617 wt(7,7) = 0.0697015 wt(7,8) = -0.00551175 wt(7,9) = 0.000856422 wt(7,10) = 0.000748415 wt(7,11) = 0.00085419 wt(7,12) = 0.000133255 wt(7,13) = -0.00024782 wt(8,1) = 0.0107957 wt(8,2) = -0.00703725 wt(8,3) = 0.00224916 wt(8,4) = -0.00214303 wt(8,5) = -0.0098033 wt(8,6) = 0.0155384 wt(8,7) = -0.00551175 wt(8,8) = 0.00425821 wt(8,9) = -0.00117158 wt(8,10) = -0.00104583 wt(8,11) = 0.000232423 wt(8,12) = -5.03151e-05 wt(8,13) = -1.04567e-05 wt(9,1) = -0.00311909 wt(9,2) = 0.00241742 wt(9,3) = -0.000777001 wt(9,4) = 0.000390261 wt(9,5) = 0.00337597 wt(9,6) = -0.00427207 wt(9,7) = 0.000856422 wt(9,8) = -0.00117158 wt(9,9) = 0.000402966 wt(9,10) = 0.000359161 wt(9,11) = -4.04778e-05 wt(9,12) = 2.81524e-05 wt(9,13) = 4.26397e-06 wt(10,1) = -0.00277854 wt(10,2) = 0.00216669 wt(10,3) = -0.000694728 wt(10,4) = 0.000344316 wt(10,5) = 0.00302563 wt(10,6) = -0.0038096 wt(10,7) = 0.000748415 wt(10,8) = -0.00104583 wt(10,9) = 0.000359161 wt(10,10) = 0.000320197 wt(10,11) = -3.62769e-05 wt(10,12) = 2.52064e-05 wt(10,13) = 4.08668e-06 wt(11,1) = 0.00113073 wt(11,2) = -0.000439147 wt(11,3) = -0.000375799 wt(11,4) = 0.000317693 wt(11,5) = -0.000621683 wt(11,6) = 0.000161267 wt(11,7) = 0.00085419 wt(11,8) = 0.000232423 wt(11,9) = -4.04778e-05 wt(11,10) = -3.62769e-05 wt(11,11) = 0.000413257 wt(11,12) = -7.32358e-06 wt(11,13) = 1.52103e-05 wt(12,1) = -0.000246653 wt(12,2) = 0.000225557 wt(12,3) = -0.000108527 wt(12,4) = 1.43079e-05 wt(12,5) = 0.000316808 wt(12,6) = -0.00045724 wt(12,7) = 0.000133255 wt(12,8) = -5.03151e-05 wt(12,9) = 2.81524e-05 wt(12,10) = 2.52064e-05 wt(12,11) = -7.32358e-06 wt(12,12) = 7.47866e-06 wt(12,13) = 1.15805e-06 wt(13,1) = 0.000106105 wt(13,2) = 0.000102069 wt(13,3) = 2.63266e-05 wt(13,4) = -7.82474e-05 wt(13,5) = 0.000142858 wt(13,6) = 0.000118812 wt(13,7) = -0.00024782 wt(13,8) = -1.04567e-05 wt(13,9) = 4.26397e-06 wt(13,10) = 4.08668e-06 wt(13,11) = 1.52103e-05 wt(13,12) = 1.15805e-06 wt(13,13) = 2.41204e-05 hdat(1) = 2.10647 hdat(2) = 0.0244452 hdat(3) = 0.076361 hdat(4) = 0.0224166 hdat(5) = 0.0343459 hdat(6) = 0.40854 hdat(7) = 0.0825449 hdat(8) = 0.708579 hdat(9) = 0.429383 hdat(10) = 0.184548 hdat(11) = 0.276277 hdat(12) = 0.000868171 hdat(13) = 0.00330593 c sum1 = 0.d0 sum2 = 0.d0 do 872 t=1,57 read(7,*) outl(t),capl(t),empl(t),hrsl(t) tshl(t) = outl(t)-theta*capl(t)-(1.d0-theta)*empl(t)-hrsl(t) sum1 = sum1 + tshl(t) sum2 = sum2 + t*tshl(t) 872 continue sum3 = 7.2055e-02*sum1-1.8797e-03*sum2 sum4 = -1.8797e-03*sum1+6.4817e-05*sum2 do 874 t=1,57 tshl(t) = tshl(t)-sum3-sum4*t 874 continue sum1 = 0.d0 sum2 = 0.d0 do 876 t=2,57 tem = tshl(t)-psi*tshl(t-1) sum1 = sum1+tem*tem sum2 = sum2+tshl(t)*tshl(t-1) 876 continue hdat(nm-1) = sum1/57.0 hdat(nm) = sum2/57.0 c army(1) = 0.00309 army(2) = 0.10663 army(3) = 0.01732 army(4) = 0.03187 mcp(1,1) = 0.962 mcp(1,2) = 0.000 mcp(1,3) = 0.019 mcp(1,4) = 0.019 c mcp(1,5) = 0.0095 mcp(2,1) = 0.000 mcp(2,2) = 0.660 mcp(2,3) = 0.000 mcp(2,4) = 0.340 c mcp(2,5) = 0.0 mcp(3,1) = 0.043 mcp(3,2) = 0.000 mcp(3,3) = 0.895 mcp(3,4) = 0.062 c mcp(3,5) = 0.0 mcp(4,1) = 0.000 mcp(4,2) = 0.115 mcp(4,3) = 0.400 mcp(4,4) = 0.485 c mcp(4,5) = 0.0 c mcp(5,1) = 0.2 c mcp(5,2) = 0.5 c mcp(5,3) = 0.0 c mcp(5,4) = 0.0 c mcp(5,5) = 0.3 c write(*,*) 'Parameters are set as follows:',param(1),param(2) write(*,*) 'Parameters are set as follows:',param(3),param(4) write(*,*) 'Parameters are set as follows:',param(5),param(6) write(*,*) 'Parameters are set as follows:',param(7),param(8) write(*,*) 'Parameters are set as follows:',param(9),param(10) write(*,*) 'Parameters are set as follows:',param(11) c xa(1) = 0.3 xa(2) = 0.4 xa(3) = 0.45 xa(4) = 0.5 xa(5) = 0.55 xa(6) = 0.6 xa(7) = 0.75 ya(1) = 0.35 ya(2) = 0.5 ya(3) = 0.65 za(1) = 0.90000E+00 za(2) = 0.100000E+01 za(3) = 0.11000E+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 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.0000001d0 maxit = 20 nx1 = nx+1 ny1 = ny+1 nz1 = nz+1 c do 40 i=1,8 do 40 j=1,8 del(i,j) =0.d0 40 continue do 50 l=1,ns ssk = 0.6 ssc = 0.12 sse = 0.40 ssh = 0.42 del(1,1) = ssk*1.0e-5 del(2,2) = ssc*1.0e-5 del(3,3) = sse*1.0e-5 del(4,4) = ssh*1.0e-5 do 60 t=1,10 d = (sse-pa)*pb dpn = 3.d0*d*(d+dsign(d,0.d0))*pb zkk = lam*ssk**theta c0 = ssc*(1.d0-ssh)**(hwgt/cwgt1) foc(1) = gamm1/gamma*sse**theta*ssc-zkk*(1.d0-ssh) foc(2) = ssc**cwgt*(1.d0-ssh)**hwgt/ome1-c0**cwgt/ome1- c dpn-gamma*ssc**cwgt1*(1.d0-ssh)**hwgt*(ssc-c0- c zkk*thet1*sse**(-theta)*ssh) foc(3) = sse*ssc+(1.d0-sse)*c0+delta*ssk+cg(l)- c zkk*ssh*sse**thet1 foc(4) = beta*(delt1+theta*zkk/ssk*sse**thet1*ssh)-1.d0 do 70 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*dssk**theta 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 = dssc**cwgt*(1.d0-dssh)**hwgt/ome1-c0**cwgt/ome1- c dpn-gamma*dssc**cwgt1*(1.d0-dssh)**hwgt*(dssc-c0- c zkk*thet1*dsse**(-theta)*dssh) foc3p = dsse*dssc+(1.d0-dsse)*c0+delta*dssk+cg(l)- c zkk*dssh*dsse**thet1 foc4p = beta*(delt1+theta*zkk/dssk*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*dssk**theta 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 = dssc**cwgt*(1.d0-dssh)**hwgt/ome1-c0**cwgt/ome1- c dpn-gamma*dssc**cwgt1*(1.d0-dssh)**hwgt*(dssc-c0- c zkk*thet1*dsse**(-theta)*dssh) foc3m = dsse*dssc+(1.d0-dsse)*c0+delta*dssk+cg(l)- c zkk*dssh*dsse**thet1 foc4m = beta*(delt1+theta*zkk/dssk*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) 70 continue call dgeco(dfoc,4,4,iwork,rcond,work) if (rcond.eq.0.0) c write(*,*) 'WARNING: singular matrix when computing ', c 'steady state values' call dgesl(dfoc,4,4,iwork,foc,0) ssk = ssk - foc(1) ssc = ssc - foc(2) sse = sse - foc(3) ssh = ssh - foc(4) 60 continue write(*,1001) ' The steady state values: ',ssk,ssc,sse,ssh ssi = delta*ssk ssn = sse ssz = 0.d0 ssf = lam*ssk**theta*sse**thet1*ssh del(1,1) = ssk*1.0e-4 del(2,2) = ssn*1.0e-4 del(3,3) = 1.0e-5 del(4,4) = 1.0e-4 del(5,5) = ssc*1.0e-4 del(6,6) = ssi*1.0e-4 del(7,7) = sse*1.0e-4 del(8,8) = ssh*1.0e-4 c d = (sse-pa)*pb pn = d*d*d+dabs(d*d*d) c0 = ssc*(1.d0-ssh)**(hwgt/cwgt1) uu = (1.d0-army(l))*((sse*ssc**cwgt*(1.d0-ssh)**hwgt+ c (1.d0-sse)*c0**cwgt)/ome1-pn) du(4) = 0.d0 do 80 i=1,8 du(i) = 0.d0 do 80 j=1,8 du2(i,j) = 0.d0 80 continue do 90 i=1,8 dssk = ssk + del(1,i) dssn = ssn + 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) d = (dsse-pa)*pb pn = d*d*d+dabs(d*d*d) c0 = (lam*dexp(dssz)*dssk**theta*dsse**thet1*dssh- c dssi-cg(l)-dsse*dssc-alpha*(dsse-dssn)*(dsse-dssn))/ c (1.d0-dsse) up = (1.d0-army(l))*((dsse*dssc**cwgt*(1.d0-dssh)**hwgt+ c (1.d0-dsse)*c0**cwgt)/ome1-pn) dssk = ssk - del(1,i) dssn = ssn - 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) d = (dsse-pa)*pb pn = d*d*d+dabs(d*d*d) c0 = (lam*dexp(dssz)*dssk**theta*dsse**thet1*dssh- c dssi-cg(l)-dsse*dssc-alpha*(dsse-dssn)*(dsse-dssn))/ c (1.d0-dsse) um = (1.d0-army(l))*((dsse*dssc**cwgt*(1.d0-dssh)**hwgt+ c (1.d0-dsse)*c0**cwgt)/ome1-pn) du(i)= 0.5*(up-um)/del(i,i) do 100 j=1,i dssk = ssk + del(1,i) + del(1,j) dssn = ssn + 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) d = (dsse-pa)*pb pn = d*d*d+dabs(d*d*d) c0 = (lam*dexp(dssz)*dssk**theta*dsse**thet1*dssh- c dssi-cg(l)-dsse*dssc-alpha*(dsse-dssn)* c (dsse-dssn))/(1.d0-dsse) upp = (1.d0-army(l))*((dsse*dssc**cwgt*(1.d0-dssh)**hwgt+ c (1.d0-dsse)*c0**cwgt)/ome1-pn) dssk = ssk + del(1,i) - del(1,j) dssn = ssn + 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) d = (dsse-pa)*pb pn = d*d*d+dabs(d*d*d) c0 = (lam*dexp(dssz)*dssk**theta*dsse**thet1*dssh- c dssi-cg(l)-dsse*dssc-alpha*(dsse-dssn)* c (dsse-dssn))/(1.d0-dsse) upm = (1.d0-army(l))*((dsse*dssc**cwgt*(1.d0-dssh)**hwgt+ c (1.d0-dsse)*c0**cwgt)/ome1-pn) dssk = ssk - del(1,i) + del(1,j) dssn = ssn - 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) d = (dsse-pa)*pb pn = d*d*d+dabs(d*d*d) c0 = (lam*dexp(dssz)*dssk**theta*dsse**thet1*dssh- c dssi-cg(l)-dsse*dssc-alpha*(dsse-dssn)* c (dsse-dssn))/(1.d0-dsse) ump = (1.d0-army(l))*((dsse*dssc**cwgt*(1.d0-dssh)**hwgt+ c (1.d0-dsse)*c0**cwgt)/ome1-pn) dssk = ssk - del(1,i) - del(1,j) dssn = ssn - 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) d = (dsse-pa)*pb pn = d*d*d+dabs(d*d*d) c0 = (lam*dexp(dssz)*dssk**theta*dsse**thet1*dssh- c dssi-cg(l)-dsse*dssc-alpha*(dsse-dssn)* c (dsse-dssn))/(1.d0-dsse) umm = (1.d0-army(l))*((dsse*dssc**cwgt*(1.d0-dssh)**hwgt+ c (1.d0-dsse)*c0**cwgt)/ome1-pn) du2(i,j) = 0.25*(upp-upm-ump+umm)/(del(i,i)*del(j,j)) du2(j,i) = du2(i,j) 100 continue 90 continue do 110 i=1,8 do 110 j=1,8 qwr(i,j) = 0.d0 110 continue do 120 j=1,8 qwr(1,j) = 0.5*du2(1,j) qwr(2,j) = 0.5*du2(2,j) qwr(3,j) = 0.5*du2(3,j) qwr(4,j) = 0.5*du2(4,j) qwr(5,j) = 0.5*du2(5,j) qwr(6,j) = 0.5*du2(6,j) qwr(7,j) = 0.5*du2(7,j) qwr(8,j) = 0.5*du2(8,j) 120 continue do 130 i=1,8 work(i) = 0.5*(du(i)-du2(i,1)*ssk-du2(i,2)*ssn-du2(i,5)* c ssc-du2(i,6)*ssi-du2(i,7)*sse-du2(i,8)*ssh) qwr(4,i) = qwr(4,i) + work(i) qwr(i,4) = qwr(i,4) + work(i) 130 continue qwr(4,4) = qwr(4,4) + uu-(work(1)+0.5*du(1))*ssk- c (work(2)+0.5*du(2))*ssn-(work(5)+0.5*du(5))*ssc- c (work(6)+0.5*du(6))*ssi-(work(7)+0.5*du(7))*sse- c (work(8)+0.5*du(8))*ssh do 140 i=1,4 do 140 j=1,4 qmat(i,j) = qwr(i,j) amat(i,j) = 0.d0 140 continue do 150 i=1,4 do 150 j=1,4 wmat(i,j) = qwr(i,j+4) bmat(i,j) = 0.d0 150 continue do 160 i=1,4 do 160 j=1,4 rmat(i,j) = qwr(i+4,j+4) 160 continue sb = dsqrt(beta) amat(1,1) = delt1*sb amat(3,3) = psi*sb amat(4,4) = sb bmat(1,2) = sb bmat(2,3) = sb c c Convert to the no-discount, no-cross product case: c do 170 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) 170 continue call dgeco(ri,ndv,ndv,iwork,rcond,work) call dgedi(ri,ndv,ndv,iwork,det,work,11) do 180 i=1,4 do 180 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) 180 continue do 190 i=1,4 do 190 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(3,j)-wmat(i,4)*rw(4,j) 190 continue do 200 i=1,4 do 200 j=1,4 amat(i,j) = amat(i,j) -bmat(i,1)*rw(1,j)-bmat(i,2)* c rw(2,j)-bmat(i,3)*rw(3,j)-bmat(i,4)*rw(4,j) 200 continue c c Doubling algorithm used to solve fixed point of Riccati equation c do 210 i=1,nsv do 210 j=1,ndv sum1 = 0.d0 do 220 k=1,ndv sum1 = sum1 + bmat(i,k)*ri(k,j) 220 continue wmat(i,j) = sum1 210 continue do 230 i=1,nsv do 230 j=1,nsv sum1 = 0.d0 do 240 k=1,ndv sum1 = sum1 + wmat(i,k)*bmat(j,k) 240 continue bb(i,j) = sum1 aa(i,j) = amat(i,j) 230 continue do 250 i=1,ndv do 250 j=1,nsv sum1 = 0.d0 do 260 k=1,nsv sum1 = sum1 + bmat(k,i)*qmat(k,j) 260 continue wmat(j,i) = sum1 250 continue do 270 i=1,ndv do 270 j=1,nsv sum1 = 0.d0 do 280 k=1,nsv sum1 = sum1 + wmat(k,i)*amat(k,j) 280 continue f0(i,j) = sum1 270 continue do 290 i=1,ndv do 290 j=1,ndv sum1 = 0.d0 do 300 k=1,nsv sum1 = sum1 + wmat(k,i)*bmat(k,j) 300 continue rb(i,j) = sum1+rmat(i,j) 290 continue call dgeco(rb,ndv,ndv,iwork,rcond,work) if (rcond.eq.0.0) c write(*,*) 'WARNING: singular matrix when computing ', c 'linear-quadratic decision function' do 310 j=1,nsv call dgesl(rb,ndv,ndv,iwork,f0(1,j),0) 310 continue d1 = 1.d0 do 320 it=1,maxit do 330 i=1,nsv do 340 j=1,nsv sum1 = 0.d0 do 350 k=1,nsv sum1 = sum1 + bb(i,k)*qmat(k,j) 350 continue ee(i,j) = sum1 340 continue ee(i,i) = ee(i,i) + 1.d0 330 continue call dgeco(ee,nsv,nsv,iwork,rcond,work) call dgedi(ee,nsv,nsv,iwork,det,work,01) do 360 i=1,nsv do 360 j=1,nsv sum1 = 0.d0 do 370 k=1,nsv sum1 = sum1 + bb(i,k)*aa(j,k) 370 continue cc(i,j) = sum1 360 continue do 380 i=1,nsv do 380 j=1,nsv sum1 = 0.d0 do 390 k=1,nsv sum1 = sum1 + aa(i,k)*ee(k,j) 390 continue dd(i,j) = sum1 380 continue do 400 i=1,nsv do 400 j=1,nsv sum1 = 0.d0 do 410 k=1,nsv sum1 = sum1 + dd(i,k)*cc(k,j) 410 continue bb(i,j) = bb(i,j)+sum1 400 continue do 420 i=1,nsv do 420 j=1,nsv sum1 = 0.d0 do 430 k=1,nsv sum1 = sum1 + ee(i,k)*aa(k,j) 430 continue cc(i,j) = sum1 420 continue do 440 i=1,nsv do 440 j=1,nsv sum1 = 0.d0 do 450 k=1,nsv sum1 = sum1 + aa(k,i)*qmat(k,j) 450 continue dd(i,j) = sum1 440 continue do 460 i=1,nsv do 460 j=1,nsv sum1 = 0.d0 do 470 k=1,nsv sum1 = sum1 + dd(i,k)*cc(k,j) 470 continue qmat(i,j) = qmat(i,j)+sum1 460 continue do 480 i=1,nsv do 480 j=1,nsv sum1 = 0.d0 do 490 k=1,nsv sum1 = sum1 + ee(i,k)*aa(k,j) 490 continue cc(i,j) = sum1 480 continue do 500 i=1,nsv do 500 j=1,nsv sum1 = 0.d0 do 510 k=1,nsv sum1 = sum1 + aa(i,k)*cc(k,j) 510 continue dd(i,j) = sum1 500 continue do 520 i=1,nsv do 520 j=1,nsv aa(i,j) = dd(i,j) 520 continue do 530 i=1,ndv do 530 j=1,nsv sum1 = 0.d0 do 540 k=1,nsv sum1 = sum1 + bmat(k,i)*qmat(k,j) 540 continue wmat(j,i) = sum1 530 continue do 550 i=1,ndv do 550 j=1,nsv sum1 = 0.d0 do 560 k=1,nsv sum1 = sum1 + wmat(k,i)*amat(k,j) 560 continue f1(i,j) = sum1 550 continue do 570 i=1,ndv do 570 j=1,ndv sum1 = 0.d0 do 580 k=1,nsv sum1 = sum1 + wmat(k,i)*bmat(k,j) 580 continue rb(i,j) = sum1+rmat(i,j) 570 continue call dgeco(rb,ndv,ndv,iwork,rcond,work) if (rcond.eq.0.0) c write(*,*) 'WARNING: singular matrix when computing ', c 'linear-quadratic decision function' do 590 j=1,nsv call dgesl(rb,ndv,ndv,iwork,f1(1,j),0) 590 continue sum1 = 0.d0 do 600 i=1,ndv do 600 j=1,nsv sum1 = sum1+(f1(i,j)-f0(i,j))*(f1(i,j)-f0(i,j)) 600 continue if (dsqrt(sum1)/float(ndv*nsv).lt.crit) go to 997 do 610 i=1,ndv do 610 j=1,nsv f0(i,j) = f1(i,j) 610 continue 320 continue 997 do 620 i=1,ndv do 620 j=1,nsv f1(i,j) = -f1(i,j)-rw(i,j) 620 continue n = 1 do 630 k=1,nz1 do 630 j=1,ny1 do 630 i=1,nx1 a0(n,l) = f1(1,1)*xa(i)+f1(1,2)*ya(j)+ c f1(1,3)*dlog(za(k))+f1(1,4) b0(n,l) = f1(3,1)*xa(i)+f1(3,2)*ya(j)+ c f1(3,3)*dlog(za(k))+f1(3,4) n = n+1 630 continue if (l.eq.3) then sski = ssk sshi = ssh ssci = ssc ssgi = ginv(l)/delta ssni = ssn endif 50 continue c do 615 i=1,npar do 625 j=1,npar+1 delp(i,j) = 0.d0 625 continue delp(i,i+1) = dmax1(dabs(param(i))*1.0e-4,1.0e-7) 615 continue do 635 it0 = 1,npar+1 alpha = param(1) + delp(1,it0) gamma = param(2) + delp(2,it0) pa = param(3) + delp(3,it0) pb = param(4) + delp(4,it0) psi = param(5) + delp(5,it0) sig = param(6) + delp(6,it0) theta = param(7) + delp(7,it0) zeta = param(8) + delp(8,it0) cg(2) = param(9) + delp(9,it0) cg(3) = param(10) + delp(10,it0) cg(4) = param(11) + delp(11,it0) sum1 = 0.d0 sum2 = 0.d0 call qgausl(ne1,-2.88d0*sig,2.88d0*sig,e1,w1) do 640 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) 640 continue do 650 i=1,ne1 wgt(i) = wgt(i)/sum1 sum2 = sum2 + wgt(i) 650 continue do 660 i=1,nx call qgausl(nxp(i),xa(i),xa(i+1),x,u) do 660 j=1,nxp(i) xx(j,i) = x(j) wx(j,i) = u(j) 660 continue do 670 i=1,ny call qgausl(nyp(i),ya(i),ya(i+1),y,v) do 670 j=1,nyp(i) yy(j,i) = y(j) wy(j,i) = v(j) 670 continue do 680 i=1,nz call qgausl(nzp(i),za(i),za(i+1),z,w) do 680 j=1,nzp(i) zz(j,i) = z(j) wz(j,i) = w(j) 680 continue ne = nx*ny*nz nxy = nx*ny sumres = 1.d0 do 690 it=1,maxit do 700 i=1,nas f(i) = 0.d0 do 700 j=1,nas df(i,j) = 0.d0 700 continue do 710 l=1,ns army1 = 1.d0-army(l) do 720 n=1,ne iz = (n-1)/nxy+1 iy = (n-(iz-1)*nxy-1)/nx+1 ix = n-(iy-1)*nx-(iz-1)*nxy n1 = nxp(ix) n2 = nyp(iy) n3 = nzp(iz) x1 = xa(ix) x2 = xa(ix+1) dx = x2-x1 y1 = ya(iy) y2 = ya(iy+1) dy = y2-y1 z1 = za(iz) z2 = za(iz+1) dz = z2-z1 nd1= ix+nx1*(iy-1)+nx1*ny1*(iz-1) nd2= nd1+1 nd4= nd1+nx1 nd3= nd4+1 nd5= nd1+nx1*ny1 nd6= nd5+1 nd8= nd5+nx1 nd7= nd8+1 do 730 i=1,n1 x(i) = xx(i,ix) u(i) = wx(i,ix) 730 continue do 740 j=1,n2 y(j) = yy(j,iy) v(j) = wy(j,iy) 740 continue do 750 k=1,n3 z(k) = zz(k,iz) w(k) = wz(k,iz) zp = (z2 - z(k))/dz zm = (z(k) - z1)/dz do 750 j=1,n2 yp = (y2 - y(j))/dy ym = (y(j) - y1)/dy do 750 i=1,n1 xp = (x2 - x(i))/dx xm = (x(i) - x1)/dx c1 = 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 n0 = b0(nd1,l)*xp*yp*zp+b0(nd2,l)*xm*yp*zp+ c b0(nd3,l)*xm*ym*zp+b0(nd4,l)*xp*ym*zp+ c b0(nd5,l)*xp*yp*zm+b0(nd6,l)*xm*yp*zm+ c b0(nd7,l)*xm*ym*zm+b0(nd8,l)*xp*ym*zm h = 1-gamm1*c1*(n0/x(i))**theta/(z(k)*lam*gamma) c0 = c1*(1.d0-h)**(hwgt/cwgt1) h1 = 1.d0-h mu = gamma*c1**cwgt1*h1**hwgt cg1 = cg(l)+zeta*dlog(z(k)) xt = delt1*x(i)+z(k)*lam*x(i)**theta*n0**thet1*h- c cg1-n0*c1-(1.d0-n0)*c0-alpha*(n0-y(j))* c (n0-y(j)) yt = n0 ii = 1 CFPP$ UNROLL do 760 m=1,nx j1 = idint( dmax1(xt-xa(m),0.d0)/(xt-xa(m))+.2 ) ii = ii + j1 760 continue ii = max(1,ii-1) jj = 1 CFPP$ UNROLL do 770 m=1,ny j1 = idint( dmax1(yt-ya(m),0.d0)/(yt-ya(m))+.2) jj = jj + j1 770 continue jj = max(1,jj-1) CFPP$ UNROLL do 780 m=1,ne1 zzt(m) = dexp(psi*dlog(z(k))+e1(m)) kk(m) = 1 CFPP$ UNROLL do 790 mm=1,nz j1 = idint( dmax1(zzt(m)-za(mm),0.d0)/ c (zzt(m)-za(mm))+.2) kk(m)= kk(m) + j1 790 continue kk(m) = max(1,kk(m)-1) 780 continue do 800 m=1,34 sum(m) = 0.d0 800 continue do 810 m=1,nas work(m) = 0.d0 wrk2(m) = 0.d0 810 continue d = (n0-pa)*pb dpn = 3.d0*d*(d+dsign(d,0.d0))*pb d2pn = 6.d0*(d+dsign(d,0.d0))*pb*pb dc0ca = -c0/c1*omega/cwgt1 dc0na = hwgt/cwgt1*c0*theta/n0 dhca = -h1/c1 dhna = -theta*h1/n0 dxca = z(k)*lam*x(i)**theta*n0**thet1*dhca-n0- c (1.d0-n0)*dc0ca dxna = z(k)*lam*x(i)**theta*n0**thet1*h*(thet1/ c n0+dhna/h)-c1+c0-(1.d0-n0)*dc0na-2.d0* c alpha*(n0-y(j)) CFPP$ UNROLL do 820 m=1,ns CFPP$ UNROLL do 820 mm=1,ne1 nct = (m-1)*ne1+mm zt = zzt(mm) x1t = xa(ii) x2t = xa(ii+1) dxt = x2t-x1t y1t = ya(jj) y2t = ya(jj+1) dyt = y2t-y1t z1t = za(kk(mm)) z2t = za(kk(mm)+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 = (kk(mm)-1)*nx1*ny1+(jj-1)*nx1+ii 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 n0t = b0(nd1t,m)*xpt*ypt*zpt+ c b0(nd2t,m)*xmt*ypt*zpt+ c b0(nd3t,m)*xmt*ymt*zpt+ c b0(nd4t,m)*xpt*ymt*zpt+ c b0(nd5t,m)*xpt*ypt*zmt+ c b0(nd6t,m)*xmt*ypt*zmt+ c b0(nd7t,m)*xmt*ymt*zmt+ c b0(nd8t,m)*xpt*ymt*zmt ht = 1.d0-gamm1*c1t*(n0t/xt)**theta/ c (gamma*zt*lam) d = (n0t-pa)*pb dpnt = 3.d0*d*(d+dsign(d,0.d0))*pb d2pnt = 6.d0*(d+dsign(d,0.d0))*pb*pb dc1tx = (-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 dc1ty = (-a0(nd1t,m)*xpt*zpt-a0(nd2t,m)*xmt*zpt c +a0(nd3t,m)*xmt*zpt+a0(nd4t,m)*xpt*zpt c -a0(nd5t,m)*xpt*zmt-a0(nd6t,m)*xmt*zmt c +a0(nd7t,m)*xmt*zmt+a0(nd8t,m)*xpt*zmt) c /dyt dn0tx = (-b0(nd1t,m)*ypt*zpt+b0(nd2t,m)*ypt*zpt c +b0(nd3t,m)*ymt*zpt-b0(nd4t,m)*ymt*zpt c -b0(nd5t,m)*ypt*zmt+b0(nd6t,m)*ypt*zmt c +b0(nd7t,m)*ymt*zmt-b0(nd8t,m)*ymt*zmt) c /dxt dn0ty = (-b0(nd1t,m)*xpt*zpt-b0(nd2t,m)*xmt*zpt c +b0(nd3t,m)*xmt*zpt+b0(nd4t,m)*xpt*zpt c -b0(nd5t,m)*xpt*zmt-b0(nd6t,m)*xmt*zmt c +b0(nd7t,m)*xmt*zmt+b0(nd8t,m)*xpt*zmt) c /dyt out = beta*zt*lam*xt**theta*n0t**thet1*ht h1t = 1.d0-ht term1 = gamma*c1t**cwgt1*h1t**hwgt*(1.d0-army(m)) term2 = beta*delt1+theta*out/xt term3 = 2.d0*alpha*beta*(n0t-yt) dc1tca = dc1tx*dxca dn0tca = dn0tx*dxca dc1tna = dc1tx*dxna+dc1ty dn0tna = dn0tx*dxna+dn0ty dhtca =-h1t*(dc1tca/c1t+theta*dn0tca/n0t-theta* c dxca/xt) dhtna =-h1t*(dc1tna/c1t+theta*dn0tna/n0t-theta* c dxna/xt) dhtcb =-h1t/c1t dhtnb =-theta*h1t/n0t dt1ca = term1*(cwgt1*dc1tca/c1t-hwgt*dhtca/h1t) dt1cb = term1*(cwgt1/c1t-hwgt*dhtcb/h1t) dt1na = term1*(cwgt1*dc1tna/c1t-hwgt*dhtna/h1t) dt1nb =-term1*hwgt*dhtnb/h1t dt2ca = theta*out/xt*(-thet1*dxca/xt+thet1* c dn0tca/n0t+dhtca/ht) dt2cb = theta*out/xt*dhtcb/ht dt2na = theta*out/xt*(-thet1*dxna/xt+thet1* c dn0tna/n0t+dhtna/ht) dt2nb = theta*out/xt*(thet1/n0t+dhtnb/ht) dt3ca = 2.d0*alpha*beta*dn0tca dt3cb = 0.d0 dt3na = 2.d0*alpha*beta*(dn0tna-1.d0) dt3nb = 2.d0*alpha*beta sum(1) = sum(1) + mcp(l,m)*wgt(mm)*term1*term2 sum(2) = sum(2) + mcp(l,m)*wgt(mm)*term1*term3 drhsc = mcp(l,m)*wgt(mm)*(dt1ca*term2+ c term1*dt2ca) drhsn = mcp(l,m)*wgt(mm)*(dt1na*term2+ c term1*dt2na) c rhs of res1,res2 c d(res1) wrt a0(nd1),b0(nd1), .... sum(3) = sum(3) + drhsc*xp*yp*zp sum(4) = sum(4) + drhsn*xp*yp*zp sum(5) = sum(5) + drhsc*xm*yp*zp sum(6) = sum(6) + drhsn*xm*yp*zp sum(7) = sum(7) + drhsc*xm*ym*zp sum(8) = sum(8) + drhsn*xm*ym*zp sum(9) = sum(9) + drhsc*xp*ym*zp sum(10)= sum(10) + drhsn*xp*ym*zp sum(11)= sum(11) + drhsc*xp*yp*zm sum(12)= sum(12) + drhsn*xp*yp*zm sum(13)= sum(13) + drhsc*xm*yp*zm sum(14)= sum(14) + drhsn*xm*yp*zm sum(15)= sum(15) + drhsc*xm*ym*zm sum(16)= sum(16) + drhsn*xm*ym*zm sum(17)= sum(17) + drhsc*xp*ym*zm sum(18)= sum(18) + drhsn*xp*ym*zm drhsc = mcp(l,m)*wgt(mm)*(dt1ca*term3+ c term1*dt3ca) drhsn = mcp(l,m)*wgt(mm)*(dt1na*term3+ c term1*dt3na) c d(res2) wrt a0(nd1),b0(nd1), .... sum(19)= sum(19) + drhsc*xp*yp*zp sum(20)= sum(20) + drhsn*xp*yp*zp sum(21)= sum(21) + drhsc*xm*yp*zp sum(22)= sum(22) + drhsn*xm*yp*zp sum(23)= sum(23) + drhsc*xm*ym*zp sum(24)= sum(24) + drhsn*xm*ym*zp sum(25)= sum(25) + drhsc*xp*ym*zp sum(26)= sum(26) + drhsn*xp*ym*zp sum(27)= sum(27) + drhsc*xp*yp*zm sum(28)= sum(28) + drhsn*xp*yp*zm sum(29)= sum(29) + drhsc*xm*yp*zm sum(30)= sum(30) + drhsn*xm*yp*zm sum(31)= sum(31) + drhsc*xm*ym*zm sum(32)= sum(32) + drhsn*xm*ym*zm sum(33)= sum(33) + drhsc*xp*ym*zm sum(34)= sum(34) + drhsn*xp*ym*zm c drhsc = mcp(l,m)*wgt(mm)*(dt1cb*term2+ c term1*dt2cb) drhsn = mcp(l,m)*wgt(mm)*(dt1nb*term2+ c term1*dt2nb) c d(res1) wrt a0(nd1t),b0(nd1t), ... wrk(nct,1) = drhsc*xpt*ypt*zpt wrk(nct,2) = drhsn*xpt*ypt*zpt wrk(nct,3) = drhsc*xmt*ypt*zpt wrk(nct,4) = drhsn*xmt*ypt*zpt wrk(nct,5) = drhsc*xmt*ymt*zpt wrk(nct,6) = drhsn*xmt*ymt*zpt wrk(nct,7) = drhsc*xpt*ymt*zpt wrk(nct,8) = drhsn*xpt*ymt*zpt wrk(nct,9) = drhsc*xpt*ypt*zmt wrk(nct,10)= drhsn*xpt*ypt*zmt wrk(nct,11)= drhsc*xmt*ypt*zmt wrk(nct,12)= drhsn*xmt*ypt*zmt wrk(nct,13)= drhsc*xmt*ymt*zmt wrk(nct,14)= drhsn*xmt*ymt*zmt wrk(nct,15)= drhsc*xpt*ymt*zmt wrk(nct,16)= drhsn*xpt*ymt*zmt drhsc = mcp(l,m)*wgt(mm)*(dt1cb*term3+ c term1*dt3cb) drhsn = mcp(l,m)*wgt(mm)*(dt1nb*term3+ c term1*dt3nb) c d(res2) wrt a0(nd1t),b0(nd1t), ... wrk(nct,17)= drhsc*xpt*ypt*zpt wrk(nct,18)= drhsn*xpt*ypt*zpt wrk(nct,19)= drhsc*xmt*ypt*zpt wrk(nct,20)= drhsn*xmt*ypt*zpt wrk(nct,21)= drhsc*xmt*ymt*zpt wrk(nct,22)= drhsn*xmt*ymt*zpt wrk(nct,23)= drhsc*xpt*ymt*zpt wrk(nct,24)= drhsn*xpt*ymt*zpt wrk(nct,25)= drhsc*xpt*ypt*zmt wrk(nct,26)= drhsn*xpt*ypt*zmt wrk(nct,27)= drhsc*xmt*ypt*zmt wrk(nct,28)= drhsn*xmt*ypt*zmt wrk(nct,29)= drhsc*xmt*ymt*zmt wrk(nct,30)= drhsn*xmt*ymt*zmt wrk(nct,31)= drhsc*xpt*ymt*zmt wrk(nct,32)= drhsn*xpt*ymt*zmt 820 continue m = (l-1)*2*na i1 = m+2*nd1-1 i2 = i1+1 i3 = m+2*nd2-1 i4 = i3+1 i5 = m+2*nd3-1 i6 = i5+1 i7 = m+2*nd4-1 i8 = i7+1 i9 = m+2*nd5-1 i10 = i9+1 i11 = m+2*nd6-1 i12 = i11+1 i13 = m+2*nd7-1 i14 = i13+1 i15 = m+2*nd8-1 i16 = i15+1 dlhsc = -mu*army1*(cwgt1/c1-hwgt*dhca/h1) dlhsn = mu*army1*hwgt*dhna/h1 work(i1) = sum(3)+dlhsc*xp*yp*zp work(i2) = sum(4)+dlhsn*xp*yp*zp work(i3) = sum(5)+dlhsc*xm*yp*zp work(i4) = sum(6)+dlhsn*xm*yp*zp work(i5) = sum(7)+dlhsc*xm*ym*zp work(i6) = sum(8)+dlhsn*xm*ym*zp work(i7) = sum(9)+dlhsc*xp*ym*zp work(i8) = sum(10)+dlhsn*xp*ym*zp work(i9) = sum(11)+dlhsc*xp*yp*zm work(i10) = sum(12)+dlhsn*xp*yp*zm work(i11) = sum(13)+dlhsc*xm*yp*zm work(i12) = sum(14)+dlhsn*xm*yp*zm work(i13) = sum(15)+dlhsc*xm*ym*zm work(i14) = sum(16)+dlhsn*xm*ym*zm work(i15) = sum(17)+dlhsc*xp*ym*zm work(i16) = sum(18)+dlhsn*xp*ym*zm tem1 = c1-c0-z(k)*lam*(x(i)/n0)**theta*thet1*h+ c 2.d0*alpha*(n0-y(j)) dlhsc =-army1*(mu-gamm1/gamma*mu*c1/h1*dhca- c gamma*c0**cwgt1*dc0ca-tem1*(cwgt1*mu/c1- c hwgt*mu/h1*dhca)-mu*(1.d0-dc0ca-z(k)* c lam*(x(i)/n0)**theta*thet1*dhca)) dlhsn =-army1*(-gamm1/gamma*mu*c1/h1*dhna-gamma* c c0**cwgt1*dc0na-d2pn+tem1*hwgt*mu/h1* c dhna-mu*(-dc0na+z(k)*lam*(x(i)/n0)** c theta*thet1*h*(theta/n0-dhna/h)+2.d0* c alpha)) wrk2(i1) = sum(19)+dlhsc*xp*yp*zp wrk2(i2) = sum(20)+dlhsn*xp*yp*zp wrk2(i3) = sum(21)+dlhsc*xm*yp*zp wrk2(i4) = sum(22)+dlhsn*xm*yp*zp wrk2(i5) = sum(23)+dlhsc*xm*ym*zp wrk2(i6) = sum(24)+dlhsn*xm*ym*zp wrk2(i7) = sum(25)+dlhsc*xp*ym*zp wrk2(i8) = sum(26)+dlhsn*xp*ym*zp wrk2(i9) = sum(27)+dlhsc*xp*yp*zm wrk2(i10) = sum(28)+dlhsn*xp*yp*zm wrk2(i11) = sum(29)+dlhsc*xm*yp*zm wrk2(i12) = sum(30)+dlhsn*xm*yp*zm wrk2(i13) = sum(31)+dlhsc*xm*ym*zm wrk2(i14) = sum(32)+dlhsn*xm*ym*zm wrk2(i15) = sum(33)+dlhsc*xp*ym*zm wrk2(i16) = sum(34)+dlhsn*xp*ym*zm CFPP$ UNROLL do 830 m=1,ns CFPP$ UNROLL do 830 mm=1,ne1 nct = (m-1)*ne1+mm j1 = (m-1)*2*na+2*ind(nct)-1 j2 = j1+1 j3 = j1+2 j4 = j3+1 j5 = j3+2*nx1 j6 = j5+1 j7 = j1+2*nx1 j8 = j7+1 j9 = j1+2*nx1*ny1 j10 = j9+1 j11 = j3+2*nx1*ny1 j12 = j11+1 j13 = j11+2*nx1 j14 = j13+1 j15 = j9+2*nx1 j16 = j15+1 work(j1) = work(j1) + wrk(nct,1) work(j2) = work(j2) + wrk(nct,2) work(j3) = work(j3) + wrk(nct,3) work(j4) = work(j4) + wrk(nct,4) work(j5) = work(j5) + wrk(nct,5) work(j6) = work(j6) + wrk(nct,6) work(j7) = work(j7) + wrk(nct,7) work(j8) = work(j8) + wrk(nct,8) work(j9) = work(j9) + wrk(nct,9) work(j10) = work(j10) + wrk(nct,10) work(j11) = work(j11) + wrk(nct,11) work(j12) = work(j12) + wrk(nct,12) work(j13) = work(j13) + wrk(nct,13) work(j14) = work(j14) + wrk(nct,14) work(j15) = work(j15) + wrk(nct,15) work(j16) = work(j16) + wrk(nct,16) wrk2(j1) = wrk2(j1) + wrk(nct,17) wrk2(j2) = wrk2(j2) + wrk(nct,18) wrk2(j3) = wrk2(j3) + wrk(nct,19) wrk2(j4) = wrk2(j4) + wrk(nct,20) wrk2(j5) = wrk2(j5) + wrk(nct,21) wrk2(j6) = wrk2(j6) + wrk(nct,22) wrk2(j7) = wrk2(j7) + wrk(nct,23) wrk2(j8) = wrk2(j8) + wrk(nct,24) wrk2(j9) = wrk2(j9) + wrk(nct,25) wrk2(j10) = wrk2(j10) + wrk(nct,26) wrk2(j11) = wrk2(j11) + wrk(nct,27) wrk2(j12) = wrk2(j12) + wrk(nct,28) wrk2(j13) = wrk2(j13) + wrk(nct,29) wrk2(j14) = wrk2(j14) + wrk(nct,30) wrk2(j15) = wrk2(j15) + wrk(nct,31) wrk2(j16) = wrk2(j16) + wrk(nct,32) 830 continue res1 = (-mu*army1+sum(1))*u(i)*v(j)*w(k) res2 = ((-c1**cwgt*h1**hwgt/ome1+c0**cwgt/ c ome1+dpn+mu*(c1-c0-z(k)*lam*(x(i)/ c n0)**theta*thet1*h+2.d0*alpha*(n0- c y(j))))*army1+sum(2))*u(i)*v(j)*w(k) f(i1) = f(i1) + res1*xp*yp*zp f(i2) = f(i2) + res2*xp*yp*zp f(i3) = f(i3) + res1*xm*yp*zp f(i4) = f(i4) + res2*xm*yp*zp f(i5) = f(i5) + res1*xm*ym*zp f(i6) = f(i6) + res2*xm*ym*zp f(i7) = f(i7) + res1*xp*ym*zp f(i8) = f(i8) + res2*xp*ym*zp f(i9) = f(i9) + res1*xp*yp*zm f(i10) = f(i10) + res2*xp*yp*zm f(i11) = f(i11) + res1*xm*yp*zm f(i12) = f(i12) + res2*xm*yp*zm f(i13) = f(i13) + res1*xm*ym*zm f(i14) = f(i14) + res2*xm*ym*zm f(i15) = f(i15) + res1*xp*ym*zm f(i16) = f(i16) + res2*xp*ym*zm c tem3 = u(i)*v(j)*w(k) do 840 mm=1,nas tem1 = tem3*work(mm) tem2 = tem3*wrk2(mm) df(i1,mm) = df(i1,mm) + tem1*xp*yp*zp df(i2,mm) = df(i2,mm) + tem2*xp*yp*zp df(i3,mm) = df(i3,mm) + tem1*xm*yp*zp df(i4,mm) = df(i4,mm) + tem2*xm*yp*zp df(i5,mm) = df(i5,mm) + tem1*xm*ym*zp df(i6,mm) = df(i6,mm) + tem2*xm*ym*zp df(i7,mm) = df(i7,mm) + tem1*xp*ym*zp df(i8,mm) = df(i8,mm) + tem2*xp*ym*zp df(i9,mm) = df(i9,mm) + tem1*xp*yp*zm df(i10,mm) = df(i10,mm) + tem2*xp*yp*zm df(i11,mm) = df(i11,mm) + tem1*xm*yp*zm df(i12,mm) = df(i12,mm) + tem2*xm*yp*zm df(i13,mm) = df(i13,mm) + tem1*xm*ym*zm df(i14,mm) = df(i14,mm) + tem2*xm*ym*zm df(i15,mm) = df(i15,mm) + tem1*xp*ym*zm df(i16,mm) = df(i16,mm) + tem2*xp*ym*zm 840 continue 750 continue 720 continue 710 continue call dgeco(df,nas,nas,iwork,rcond,work) if (rcond.eq.0.0) c write(*,*) 'WARNING: singular matrix when solving df*x=f' call dgesl(df,nas,nas,iwork,f,0) sum3 = 0.d0 do 850 j=1,ns do 850 i=1,na i1 = (j-1)*na*2+2*i-1 tem1 = f(i1) tem2 = f(i1+1) a1(i,j) = a0(i,j) - tem1 b1(i,j) = b0(i,j) - tem2 sum3 = sum3 + tem1*tem1+tem2*tem2 850 continue sumres = dsqrt(sum3)/float(nas) write(*,*) ' at iteration ', it, ' the residual is ', sumres do 860 j=1,ns do 860 i=1,na a0(i,j) = a1(i,j) b0(i,j) = b1(i,j) 860 continue if (sumres.lt.crit) go to 999 690 continue 999 continue c c simulate model time series c cap = sski capg = ssgi elag = ssni tech = 1.d0 olag = lam*cap**theta*elag**thet1*sshi c1lag = ssci tlag = 1.d0 c0lag = c1lag*(1.d0-sshi)**(hwgt/cwgt1) sum1 = 0.d0 sum2 = 0.d0 sum3 = 0.d0 sum4 = 0.d0 sum5 = 0.d0 sum6 = 0.d0 sum7 = 0.d0 sum8 = 0.d0 sum9 = 0.d0 sum10 = 0.d0 sum11 = 0.d0 sum12 = 0.d0 sum13 = 0.d0 sum14 = 0.d0 sum15 = 0.d0 sum16 = 0.d0 sum17 = 0.d0 sum18 = 0.d0 do 875 i=1,ns st(i) = 0.d0 875 continue do 880 t=1,nt l = istat(t) st(l) = 1.d0 cg1 = cg(l)+zeta*dlog(tech) c c find cap,elag,tech on grid c ii = 1 do 890 m=1,nx j1 = idint( dmax1(cap-xa(m),0.d0)/(cap-xa(m))+.2 ) ii = ii + j1 890 continue ii = max(1,ii-1) jj = 1 do 900 m=1,ny j1 = idint( dmax1(elag-ya(m),0.d0)/(elag-ya(m))+.2) jj = jj + j1 900 continue jj = max(1,jj-1) k1 = 1 do 910 m=1,nz j1 = idint( dmax1(tech-za(m),0.d0)/(tech-za(m))+.2) k1 = k1 + j1 910 continue k1 = max(1,k1-1) xp = (xa(ii+1)-cap)/(xa(ii+1)-xa(ii)) xm = (cap-xa(ii))/(xa(ii+1)-xa(ii)) yp = (ya(jj+1)-elag)/(ya(jj+1)-ya(jj)) ym = (elag-ya(jj))/(ya(jj+1)-ya(jj)) zp = (za(k1+1)-tech)/(za(k1+1)-za(k1)) zm = (tech-za(k1))/(za(k1+1)-za(k1)) nd1 = (k1-1)*nx1*ny1+(jj-1)*nx1+ii nd2 = nd1+1 nd4 = nd1+nx1 nd3 = nd4+1 nd5 = nd1+nx1*ny1 nd6 = nd5+1 nd8 = nd5+nx1 nd7 = nd8+1 c1 = 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 n0 = b0(nd1,l)*xp*yp*zp+b0(nd2,l)*xm*yp*zp+ c b0(nd3,l)*xm*ym*zp+b0(nd4,l)*xp*ym*zp+ c b0(nd5,l)*xp*yp*zm+b0(nd6,l)*xm*yp*zm+ c b0(nd7,l)*xm*ym*zm+b0(nd8,l)*xp*ym*zm h = 1.d0-gamm1*c1*(n0/cap)**theta/(gamma*tech*lam) c0 = c1*(1.d0-h)**(hwgt/cwgt1) capp = cap-capg out = tech*lam*cap**theta*n0**thet1*h c c calculate moments c sum1 = sum1 + capp/out sum2 = sum2 + cg1/out*st(2) sum3 = sum3 + cg1/out*st(3) sum4 = sum4 + cg1/out*st(4) sum5 = sum5 + n0*st(2) sum6 = sum6 + n0*st(3) sum7 = sum7 + n0*st(4) tem = (c1*n0+(1.d0-n0)*c0)/out sum8 = sum8 + tem sum9 = sum9 + h sum10 = sum10 + h*h sum11 = sum11 + n0*n0 tem = dlog(tech)-psi*dlog(tlag) sum12 = sum12 + tem*tem sum13 = sum13 + dlog(tech)*dlog(tlag) c c update states c cap = out+delt1*cap-n0*c1-(1.d0-n0)*c0-cg1- c alpha*(n0-elag)*(n0-elag) elag = n0 c1lag= c1 c0lag= c0 olag = out tlag = tech tech = dexp(psi*dlog(tech)+sig*eps(t)) capg = delt1*capg+ginv(l) st(l)= 0.d0 880 continue hmod(1) = sum1/float(nt) hmod(2) = sum2/float(nt) hmod(3) = sum3/float(nt) hmod(4) = sum4/float(nt) hmod(5) = sum5/float(nt) hmod(6) = sum6/float(nt) hmod(7) = sum7/float(nt) hmod(8) = sum8/float(nt) hmod(9) = sum9/float(nt) hmod(10)= sum10/float(nt) hmod(11)= sum11/float(nt) hmod(12)= sum12/float(nt) hmod(13)= sum13/float(nt) obj = 0.d0 do 920 i=1,nm sum1 = 0.d0 do 930 k=1,nm sum1 = sum1 + wt(i,k)*(hdat(k)-hmod(k)) 930 continue obj = obj + (hdat(i)-hmod(i))*sum1 920 continue write(*,1002) ' Objective function is: ',obj write(*,*) ' Moments in the data and model ' do 940 i=1,nm write(*,*) hdat(i),hmod(i) 940 continue if (it0.eq.1) objhld = obj if (it0.gt.1) grad(it0-1) = (obj-objhld)/delp(it0-1,it0) if ((it0.eq.1).and.(isterr.eq.1)) then do 960 i=1,nm hmhld(i) = hmod(i) 960 continue endif if ((it0.gt.1).and.(isterr.eq.1)) then do 970 i=1,nm grse(i,it0-1) = (hmod(i)-hmhld(i))/delp(it0-1,it0) 970 continue endif 635 continue if (isterr.eq.1) then do 980 i=1,nm do 980 j=1,npar sum1 = 0.d0 do 990 k=1,nm sum1 = sum1 + wt(i,k)*grse(k,j) 990 continue ff(i,j) = sum1 980 continue do 1010 i=1,npar do 1010 j=1,npar sum1 = 0.d0 do 1020 k=1,nm sum1 = sum1 + grse(k,i)*ff(k,j) 1020 continue gg(i,j) = 0.016708*sum1 1010 continue do 1030 i=1,npar se(i) = dsqrt(gg(i,i)) 1030 continue endif c c var = 1/57*grse'*(20/21)*wt*grse (nm x npar)'*nm x nm *(nm x npar) c se = sqrt(diag(var)) npar x 1 c obj = objhld write(*,*) 'Parameters and gradients (and standard errors): ' do 950 i=1,npar if (isterr.eq.0) write(*,*) param(i),grad(i) if (isterr.eq.1) write(*,*) param(i),grad(i),se(i) 950 continue close(5) close(7) 1001 format (A27,4(F8.3,1X)) 1002 format (A27,E10.5,1X,2(F8.3,1X)) return end