program intang ! ! DESCRIPTION: This program computes the equilibrium for the stochastic ! version of the model in McGrattan and Prescott, Staff ! Report 369. See notes.tex and comments below for details. ! ! ! MODEL Households solve: ! ECONOMY: ! max E sum beta^t pi[t] U(c[t],1-h[t]) N[t] ! {c,xb,xi ! kbf,kbi} ! hbf,hbi} ! ! subject to ! ! c[t]+xb[t]+q[t]xi[t] ! <= sum rb[t]kb[t]+ri[t]ki[t]+w[t]h[t] ! + nonbusiness income less investment ! -all taxes + transfers ! ! kb[t+1] = ((1-deltb) kb[t] + xb[t]) / (1+eta) ! ki[t+1] = ((1-delti) ki[t] + xi[t]) / (1+eta) ! ! Rental and Wage rates: marginal products from ! ! yb = F(kbf, ki,Z1 hbf), Z1 = z1*(1+gamma)^t ! xi = G(kbi, ki,Z2 hbi), Z2 = z2*(1+gamma)^t ! h = hbf+hbi+hn ! kb = kbf+kbi ! ! Resource Constraint: ! ! c + xb + q*xi + g = y + (yn-xn) ! ! Processes for exogenous states: ! ! i is the index for a Markov chain with transition pi ! i is in the set {1,...,ns} ! ! Exogenous variables (functions of i): ! ! technology shocks (z1,z2) ! fiscal shocks (g, tax rates) ! nonbusiness activity (yn, xn, hn) ! ! ! SOLUTION Let x = tangible business capital ! y = intangible business capital ! i = index of exogenous states. ! ! Represent kbf and hbf as weighted sums of bilinear ! basis functions: ! ! kbf(x,y,i) = sum_j alpk_j^i N_j(x,y) ! hbf(x,y,i) = sum_j alph_j^i N_j(x,y) ! ! where ! N_j(x,y) = (x -x_{a-1}) (y -y_{b-1}) ! -------------* ------------- ! (x_a-x_{a-1}) (y_b-y_{b-1}) ! ! on [x_{a-1},x_a] x [y_{b-1},y_b] ! ! = (x_{a+1}-x ) (y -y_{b-1}) ! -------------* ------------- ! (x_{a+1}-x_a) (y_b-y_{b-1}) ! ! on [x_a,x_{a+1}] x [y_{b-1},y_b] ! ! = (x -x_{a-1}) (y_{b+1}-y ) ! -------------* ------------- ! (x_a-x_{a-1}) (y_{b+1}-y_b) ! ! on [x_{a-1},x_a] x [y_b,y_{b+1}] ! ! = (x_{a+1}-x ) (y_{b+1}-y ) ! -------------* ------------- ! (x_{a+1}-x_a) (y_{b+1}-y_b) ! ! on [x_a,x_{a+1}] x [y_b,y_{b+1}] ! ! = 0, elsewhere ! ! where node j is at the center of the 4 rectangles in ! which N_j(x,y) is nonzero, i.e.: ! ! --*-------------*-------------*-- ! y_{b+1}| | | ! | | | ! | | | ! | | | ! --*------------j*-------------*-- ! y_b| | | ! | | | ! | | | ! | | | ! --*-------------*-------------* ! y_{b-1},x_{a-1} |x_a |x_{a+1} ! ! * = node ! *j= center node ! ! Plug approximate kbf(x,y,i) and hbf(x,y,i) ! into the intertemporal equations and set them equal to ! 0. Assume that these equations can be represented as ! ! R(x,y,i;alp) = 0. ! ! Then the approximate solutions computed here are ! the functions given above with the alp's chosen to ! satisfy the integral equations: ! ! int R(x,y,i;alp) N_j(x,y) dx dy = 0 ! x,y ! ! for all i,j ! ! FUNCTIONAL U(c,1-h) = log(c) + psi log(1-h) ! FORMS: F(kbf,ki,hbf) = kbf^thet1 *ki^phi1 * hbf^(1-thet1-phi1) ! G(kbi,ki,hbi) = kbi^thet2 *ki^phi2 * hbi^(1-thet2-phi2) ! ! INPUT FILE: intang.inp assumes the following order for inputs: ! ! beta : discount factor ! deltb : depreciation of business tangible capital ! delti : depreciation of business intangible capital ! eta : growth of population ! gamma : growth of technology ! thet1 : share of business tangible capital, F ! phi1 : share of business intangible capital, F ! thet2 : share of business tangible capital, G ! phi2 : share of business intangible capital, G ! psi : preference parameter, weight on leisure ! chi : fraction of intangible shareholders finance ! ns : number of exogenous states ! z1 : level of technology, F (1 x ns) ! z2 : level of technology, G (1 x ns) ! tauc : tax rate on consumption (1 x ns) ! tauh : tax rate on hours (1 x ns) ! taux : tax rate on investment (1 x ns) ! tauk : tax rate on property (1 x ns) ! taup : tax rate on profits (1 x ns) ! taud : tax rate on dividends (1 x ns) ! g : government spending (1 x ns) ! yn : nonbusiness value added (1 x ns) ! xn : nonbusiness investment (1 x ns) ! hn : nonbusiness hours (1 x ns) ! pi' : transpose transition matrix (ns x ns) ! (columns sum to 1) ! x : grid on business tangible capital (nx x 1) ! y : grid on business intangible capital (ny x 1) ! kbf : guess for kbf (nx*ny*ns x 1) ! hbf : guess for hbf (nx*ny*ns x 1) ! icomp : 1=compute equilibrium, 0 otherwise ! isimul : 1=simulate time series, 0 otherwise ! nsim : number of simulations ! nt : number of time periods per simulation ! idraw : 1=draw random realizations, 0 otherwise ! iseed : integer seed for random number generator ! x0,y0 : if idraw=0, initial state and indices for ! ind : exogenous states (1 x nt) ! ! Ellen McGrattan, 8-20-07 ! Revised, ERM, 10-20-07 implicit none integer, parameter :: nx=4,ny=3,nxy=12,ns=26,na=624,nq=3, & msim=100,mt=2000 real, parameter :: crit= 0.0000000001 real, dimension (nx-1,nq) :: ax,wx real, dimension (ny-1,nq) :: ay,wy real, dimension (nx) :: xa real, dimension (ny) :: ya real, dimension (ns) :: z1,z2,tauc,tauh,taux,tauk,taup,taud,g, & yn,xn,hn real, dimension (ns,ns) :: pi = 0.0 real, dimension (nxy,ns) :: kbf0,kbf1,hbf0,hbf1 real, dimension (na,1) :: r,wrk1,wrk2,sol real, dimension (na,na) :: dr real, dimension (msim) :: x0,y0 real :: beta,beth,deltb,delti,eta,gamma,grate,psi, & thet1,phi1,lshare1,thet2,phi2,lshare2,chi, & kb,kbt,x,xt,wgtx,x1,x2,dx,x1t,x2t,dxt, & ki,kit,y,yt,wgty,y1,y2,dy,y1t,y2t,dyt, & kbf,hbf,kbft,hbft,kbi,yb,qxi,hbi,xi, & q,h,c,mu,xb,d,kbit,ybt,qxit,hbit,xit, & qt,ht,ct,rbt,rit,mut,tdx,tdp,tdph,tdxt, & tdpt,tdpht,basis1,basis2,basis3,basis4, & basis1t,basis2t,basis3t,basis4t, & tkbftxt,tkbftyt,thbftxt,thbftyt, & tkbftkbf,thbftkbf,tkbfthbf,thbfthbf, & tybtkbft,tqxitkbft,thbitkbft,txitkbft, & tqtkbft,thtkbft,tctkbft,tybthbft, & tqxithbft,thbithbft,txithbft,tqthbft, & ththbft,tcthbft,tybtkbf,tqxitkbf,thbitkbf, & txitkbf,tqtkbf,thtkbf,tctkbf,tkbitkbf, & tybthbf,tqxithbf,thbithbf,txithbf,tqthbf, & ththbf,tcthbf,tkbithbf,tybkbf,tqxikbf, & thbikbf,txikbf,tqkbf,thkbf,tckbf,txbkbf, & txtkbf,tytkbf,tybhbf,tqxihbf,thbihbf, & txihbf,tqhbf,thhbf,tchbf,txbhbf,txthbf, & tythbf,tmutkbft,trbtkbft,tritkbft, & tmuthbft,trbthbft,trithbft,tmutkbf, & trbtkbf,tritkbf,tmuthbf,trbthbf,trithbf, & tmukbf,tmuhbf,sum1,sum2,tsum1k,tsum2k, & tsum1h,tsum2h,tsum1kt,tsum2kt,tsum1ht, & tsum2ht,res1,res2,dres1k,dres2k,dres1h, & dres2h,eps,useed,ustart,rn,uni integer, dimension (na,1) :: iwrk integer, dimension (na+1,1) :: ipr integer, dimension (msim,mt):: ind integer :: init,info,ix,iy,ixt,iyt,i,j,k,l,m, & n,ne,icomp,isimul,idraw,nsim,nt,iseed, & nd1,nd2,nd3,nd4,nd1t,nd2t,nd3t,nd4t, & l1,l2,l3,l4,l5,l6,l7,l8, & m1,m2,m3,m4,m5,m6,m7,m8 ! ! ! Parameters: ! open(unit=5, file='intang.inp') open(unit=7, file='intang.nxt') open(unit=8, file='intang.dat') open(unit=9, file='ifunc.dat') read(5,*) beta read(5,*) deltb read(5,*) delti read(5,*) eta read(5,*) gamma read(5,*) thet1 read(5,*) phi1 read(5,*) thet2 read(5,*) phi2 read(5,*) psi read(5,*) chi read(5,*) i if (i /= ns) then write(*,*) 'ERROR: wrong number of exogenous states. Edit intang.f90' write(*,*) ' and change parameter ns.' stop endif read(5,*) z1,z2,tauc,tauh,taux,tauk,taup,taud,g,yn,xn,hn read(5,*) pi read(5,*) xa read(5,*) ya do i=1,nxy read(5,*) (kbf0(i,j),j=1,ns) end do do i=1,nxy read(5,*) (hbf0(i,j),j=1,ns) end do read(5,*) icomp read(5,*) isimul if (isimul == 1) then read(5,*) nsim if (nsim > msim) then write(*,*) 'ERROR: Parameter governing maximum simulations is too ' write(*,*) ' small. Edit file and increase parameter msim.' stop endif read(5,*) nt if (nt > mt) then write(*,*) 'ERROR: Parameter governing maximum simulation periods is' write(*,*) ' too small. Edit file and increase parameter mt.' stop endif read(5,*) idraw read(5,*) iseed if (idraw == 0) then do i=1,nsim read(5,*) x0(i),y0(i),(ind(i,j),j=1,nt) end do endif endif ! ! Intermediate parameters: ! beth = beta/(1.0+gamma) grate = (1.0+gamma)*(1.0+eta) lshare1 = 1.0-thet1-phi1 lshare2 = 1.0-thet2-phi2 ! ! Quadarature abscissas and weights: ! do i=1,nx-1 call qgausl (nq,xa(i),xa(i+1),ax(i,:),wx(i,:)) end do do i=1,ny-1 call qgausl (nq,ya(i),ya(i+1),ay(i,:),wy(i,:)) end do ! ! Compute equilibrium ! if (icomp == 1) then ! ! Solve the fixed point problem: int R(x,y,i;alp) dx dy = 0 ! for alp using a Newton method. ! write(*,*) 'Computing decision functions .... ' ne = (nx-1)*(ny-1) newton: do ! ! Derive residuals (r) of the first order conditions and their ! derivatives (dr) with respect to the unknowns in alpha. ! r = 0.0 dr = 0.0 do l=1,ns tdx = (1.0-taud(l))*(1.0+taux(l)) tdp = (1.0-taud(l))*(1.0-taup(l)) tdph = chi*(1.0-taud(l))*(1.0-taup(l))+(1.0-chi)*(1.0-tauh(l)) do n=1,ne ! ! Element n is [xa(ix),xa(ix+1)] x [ya(iy),ya(iy+1)] ! ix = mod(n-1,nx-1)+1 iy = (mod(n-1,(nx-1)*(ny-1))/(nx-1))+1 x1 = xa(ix) x2 = xa(ix+1) dx = x2-x1 y1 = ya(iy) y2 = ya(iy+1) dy = y2-y1 nd1 = ix+nx*(iy-1) nd2 = nd1+1 nd4 = nd1+nx nd3 = nd4+1 ! ! Compute r at all quadrature points on the element. ! do j=1,nq y = ay(iy,j) wgty = wy(iy,j) do i=1,nq x = ax(ix,i) wgtx = wx(ix,i) basis1 = (x2-x)*(y2-y)/(dx*dy) basis2 = (x-x1)*(y2-y)/(dx*dy) basis3 = (x-x1)*(y-y1)/(dx*dy) basis4 = (x2-x)*(y-y1)/(dx*dy) ! ! Compute decisions using finite element approximation. ! kbf = kbf0(nd1,l)*basis1 + kbf0(nd2,l)*basis2 & +kbf0(nd3,l)*basis3 + kbf0(nd4,l)*basis4 hbf = hbf0(nd1,l)*basis1 + hbf0(nd2,l)*basis2 & +hbf0(nd3,l)*basis3 + hbf0(nd4,l)*basis4 ! ! Back out other decisions using static equations ! kb = x ki = y kbi = kb-kbf yb = kbf**thet1*ki**phi1*(z1(l)*hbf)**lshare1 qxi = thet1*yb*kbi/(thet2*kbf) hbi = lshare2*qxi*hbf/(lshare1*yb) xi = kbi**thet2*ki**phi2*(z2(l)*hbi)**lshare2 q = qxi/xi h = hbf+hbi+hn(l) c = (1.0-tauh(l))*lshare1*yb*(1.0-h)/ & (psi*(1.0+tauc(l))*hbf) xb = yb+yn(l)-xn(l)-g(l)-c mu = 1.0/((1+tauc(l))*c) kbt = ((1.0-deltb)*kb+xb)/grate kit = ((1.0-delti)*ki+xi)/grate ! ! Derivatives to be used later ! tybkbf = thet1*yb/kbf tqxikbf = qxi*(tybkbf/yb-1.0/kbi-1.0/kbf) thbikbf = hbi*(tqxikbf/qxi-tybkbf/yb) txikbf = xi*(-thet2/kbi+lshare2*thbikbf/hbi) tqkbf = (tqxikbf-q*txikbf)/xi thkbf = thbikbf tckbf = c*(tybkbf/yb-thkbf/(1.0-h)) txbkbf = tybkbf-tckbf txtkbf = txbkbf/grate tytkbf = txikbf/grate tmukbf = -mu*tckbf/c tybhbf = lshare1*yb/hbf tqxihbf = qxi*tybhbf/yb thbihbf = hbi*(tqxihbf/qxi+1.0/hbf-tybhbf/yb) txihbf = lshare2*xi*thbihbf/hbi tqhbf = (tqxihbf-q*txihbf)/xi thhbf = 1.0+thbihbf tchbf = c*(tybhbf/yb-thhbf/(1.0-h)-1.0/hbf) txbhbf = tybhbf-tchbf txthbf = txbhbf/grate tythbf = txihbf/grate tmuhbf = -mu*tchbf/c ! ! Update states and find tomorrow's state on the grid. ! xt = kbt yt = kit ixt = 1 do m = 1,nx-1 if (xt > xa(m)) ixt = m end do x1t = xa(ixt) x2t = xa(ixt+1) dxt = x2t-x1t iyt = 1 do m = 1,ny-1 if (yt > ya(m)) iyt = m end do y1t = ya(iyt) y2t = ya(iyt+1) dyt = y2t-y1t nd1t = ixt+nx*(iyt-1) nd2t = nd1t+1 nd4t = nd1t+nx nd3t = nd4t+1 basis1t = (x2t-xt)*(y2t-yt)/(dxt*dyt) basis2t = (xt-x1t)*(y2t-yt)/(dxt*dyt) basis3t = (xt-x1t)*(yt-y1t)/(dxt*dyt) basis4t = (x2t-xt)*(yt-y1t)/(dxt*dyt) ! ! Initialize work vectors and sums. ! wrk1 = 0.0 wrk2 = 0.0 sum1 = 0.0 sum2 = 0.0 tsum1k = 0.0 tsum2k = 0.0 tsum1h = 0.0 tsum2h = 0.0 ! ! Compute the following 2 sums: ! betahat sum_s' pi(s'|s){mu'*Rb} ! betahat sum_s' pi(s'|s){mu'*Ru} ! do m=1,ns tdxt = (1.0-taud(m))*(1.0+taux(m)) tdpt = (1.0-taud(m))*(1.0-taup(m)) tdpht = chi*(1.0-taud(m))*(1.0-taup(m))+ & (1.0-chi)*(1.0-tauh(m)) ! ! Next period decisions using finite element approximation. ! kbft = kbf0(nd1t,m)*basis1t + kbf0(nd2t,m)*basis2t+ & kbf0(nd3t,m)*basis3t + kbf0(nd4t,m)*basis4t hbft = hbf0(nd1t,m)*basis1t + hbf0(nd2t,m)*basis2t+ & hbf0(nd3t,m)*basis3t + hbf0(nd4t,m)*basis4t ! ! Derivatives of variables with respect to coefficients ! in finite element approximation needed. ! tkbftxt = -(kbf0(nd1t,m)*basis1t+ & kbf0(nd4t,m)*basis4t)/(x2t-xt)+ & (kbf0(nd2t,m)*basis2t+ & kbf0(nd3t,m)*basis3t)/(xt-x1t) tkbftyt = -(kbf0(nd1t,m)*basis1t+ & kbf0(nd2t,m)*basis2t)/(y2t-yt)+ & (kbf0(nd3t,m)*basis3t+ & kbf0(nd4t,m)*basis4t)/(yt-y1t) thbftxt = -(hbf0(nd1t,m)*basis1t+ & hbf0(nd4t,m)*basis4t)/(x2t-xt)+ & (hbf0(nd2t,m)*basis2t+ & hbf0(nd3t,m)*basis3t)/(xt-x1t) thbftyt = -(hbf0(nd1t,m)*basis1t+ & hbf0(nd2t,m)*basis2t)/(y2t-yt)+ & (hbf0(nd3t,m)*basis3t+ & hbf0(nd4t,m)*basis4t)/(yt-y1t) tkbftkbf = tkbftxt*txtkbf+tkbftyt*tytkbf thbftkbf = thbftxt*txtkbf+thbftyt*tytkbf tkbfthbf = tkbftxt*txthbf+tkbftyt*tythbf thbfthbf = thbftxt*txthbf+thbftyt*tythbf ! ! Back out other variables. ! kbit = kbt-kbft ybt = kbft**thet1*kit**phi1*(z1(m)*hbft)**lshare1 qxit = thet1*ybt*kbit/(thet2*kbft) hbit = lshare2*qxit*hbft/(lshare1*ybt) xit = kbit**thet2*kit**phi2*(z2(m)*hbit)**lshare2 qt = qxit/xit ht = hbft+hbit+hn(m) ct = (1.0-tauh(m))*lshare1*ybt*(1.0-ht)/ & (psi*(1.0+tauc(m))*hbft) rbt = (tdpt*(thet1*ybt/kbft-tauk(m))+(1.0-deltb)*tdxt+ & deltb*(1.0-taud(m))*taup(m))/tdx rit = (tdpt*(phi1*ybt+phi2*qxit)/kit+ & (1.0-delti)*tdpht*qt)/(q*tdph) mut = 1.0/((1+tauc(m))*ct) sum1 = sum1 +beth*pi(l,m)*mut*rbt/mu sum2 = sum2 +beth*pi(l,m)*mut*rit/mu ! ! Compute derivatives of sum1 with respect to the ! coefficients in the finite element approximation. ! tybtkbft = thet1*ybt/kbft tqxitkbft = qxit*(tybtkbft/ybt-1.0/kbit-1.0/kbft) thbitkbft = hbit*(tqxitkbft/qxit-tybtkbft/ybt) txitkbft = xit*(-thet2/kbit+lshare2*thbitkbft/hbit) tqtkbft = (tqxitkbft-qt*txitkbft)/xit thtkbft = thbitkbft tctkbft = ct*(tybtkbft/ybt-thtkbft/(1.0-ht)) tkbitkbf = txtkbf-tkbftkbf tybtkbf = ybt*(thet1*tkbftkbf/kbft+phi1*tytkbf/kit+ & lshare1*thbftkbf/hbft) tqxitkbf = qxit*(tybtkbf/ybt+tkbitkbf/kbit-tkbftkbf/kbft) thbitkbf = hbit*(tqxitkbf/qxit+thbftkbf/hbft-tybtkbf/ybt) txitkbf = xit*(thet2*tkbitkbf/kbit+phi2*tytkbf/kit+ & lshare2*thbitkbf/hbit) tqtkbf = (tqxitkbf-qt*txitkbf)/xit thtkbf = thbftkbf+thbitkbf tctkbf = ct*(tybtkbf/ybt-thtkbf/(1.0-ht)-thbftkbf/hbft) tybthbft = lshare1*ybt/hbft tqxithbft = qxit*tybthbft/ybt thbithbft = hbit*(tqxithbft/qxit+1.0/hbft-tybthbft/ybt) txithbft = xit*(lshare2*thbithbft/hbit) tqthbft = (tqxithbft-qt*txithbft)/xit ththbft = 1.0+thbithbft tcthbft = ct*(tybthbft/ybt-ththbft/(1.0-ht)-1.0/hbft) tkbithbf = txthbf-tkbfthbf tybthbf = ybt*(thet1*tkbfthbf/kbft+phi1*tythbf/kit+ & lshare1*thbfthbf/hbft) tqxithbf = qxit*(tybthbf/ybt+tkbithbf/kbit-tkbfthbf/kbft) thbithbf = hbit*(tqxithbf/qxit+thbfthbf/hbft-tybthbf/ybt) txithbf = xit*(thet2*tkbithbf/kbit+phi2*tythbf/kit+ & lshare2*thbithbf/hbit) tqthbf = (tqxithbf-qt*txithbf)/xit ththbf = thbfthbf+thbithbf tcthbf = ct*(tybthbf/ybt-ththbf/(1.0-ht)-thbfthbf/hbft) tmutkbft = -mut*tctkbft/ct tmutkbf = -mut*tctkbf/ct trbtkbft = tdpt*thet1*(tybtkbft-ybt/kbft)/(tdx*kbft) trbtkbf = tdpt*thet1*(tybtkbf-ybt*tkbftkbf/kbft)/(tdx*kbft) tritkbft = (tdpt*(phi1*tybtkbft+phi2*tqxitkbft)/kit+ & (1.0-delti)*tdpht*tqtkbft)/(q*tdph) tritkbf = tdpt*(phi1*tybtkbf+phi2*tqxitkbf)/ & (tdph*q*kit)-tdpt*(phi1*ybt+phi2*qxit)* & (tqkbf/q+tytkbf/kit)/(tdph*q*kit)+ & (1.0-delti)*tdpht*(tqtkbf-tqkbf*qt/q)/(q*tdph) tmuthbft = -mut*tcthbft/ct tmuthbf = -mut*tcthbf/ct trbthbft = tdpt*thet1*tybthbft/(tdx*kbft) trbthbf = tdpt*thet1*(tybthbf-ybt*tkbfthbf/kbft)/(tdx*kbft) trithbft = (tdpt*(phi1*tybthbft+phi2*tqxithbft)/kit+ & (1.0-delti)*tdpht*tqthbft)/(q*tdph) trithbf = tdpt*(phi1*tybthbf+phi2*tqxithbf)/ & (tdph*q*kit)-tdpt*(phi1*ybt+phi2*qxit)* & (tqhbf/q+tythbf/kit)/(tdph*q*kit)+ & (1.0-delti)*tdpht*(tqthbf-tqhbf*qt/q)/(q*tdph) tsum1k = tsum1k +beth*pi(l,m)*(tmutkbf*rbt+mut*trbtkbf- & rbt*mut*tmukbf/mu)/mu tsum2k = tsum2k +beth*pi(l,m)*(tmutkbf*rit+mut*tritkbf- & rit*mut*tmukbf/mu)/mu tsum1h = tsum1h +beth*pi(l,m)*(tmuthbf*rbt+mut*trbthbf- & rbt*mut*tmuhbf/mu)/mu tsum2h = tsum2h +beth*pi(l,m)*(tmuthbf*rit+mut*trithbf- & rit*mut*tmuhbf/mu)/mu tsum1kt = beth*pi(l,m)*(tmutkbft*rbt+mut*trbtkbft)/mu tsum2kt = beth*pi(l,m)*(tmutkbft*rit+mut*tritkbft)/mu tsum1ht = beth*pi(l,m)*(tmuthbft*rbt+mut*trbthbft)/mu tsum2ht = beth*pi(l,m)*(tmuthbft*rit+mut*trithbft)/mu ! ! Hold results in work vectors for later. ! m1 = 2*nd1t-1+(m-1)*2*nxy m2 = m1+1 m3 = m1+2 m4 = m3+1 m5 = m3+2*nx m6 = m5+1 m7 = m1+2*nx m8 = m7+1 wrk1(m1,1)= wrk1(m1,1) + tsum1kt * basis1t wrk1(m2,1)= wrk1(m2,1) + tsum1ht * basis1t wrk1(m3,1)= wrk1(m3,1) + tsum1kt * basis2t wrk1(m4,1)= wrk1(m4,1) + tsum1ht * basis2t wrk1(m5,1)= wrk1(m5,1) + tsum1kt * basis3t wrk1(m6,1)= wrk1(m6,1) + tsum1ht * basis3t wrk1(m7,1)= wrk1(m7,1) + tsum1kt * basis4t wrk1(m8,1)= wrk1(m8,1) + tsum1ht * basis4t wrk2(m1,1)= wrk2(m1,1) + tsum2kt * basis1t wrk2(m2,1)= wrk2(m2,1) + tsum2ht * basis1t wrk2(m3,1)= wrk2(m3,1) + tsum2kt * basis2t wrk2(m4,1)= wrk2(m4,1) + tsum2ht * basis2t wrk2(m5,1)= wrk2(m5,1) + tsum2kt * basis3t wrk2(m6,1)= wrk2(m6,1) + tsum2ht * basis3t wrk2(m7,1)= wrk2(m7,1) + tsum2kt * basis4t wrk2(m8,1)= wrk2(m8,1) + tsum2ht * basis4t end do ! ! The wrk vectors contain derivatives of the residuals. ! l1 = 2*nd1-1+(l-1)*2*nxy l2 = l1+1 l3 = l1+2 l4 = l3+1 l5 = l3+2*nx l6 = l5+1 l7 = l1+2*nx l8 = l7+1 dres1k = tsum1k dres1h = tsum1h dres2k = tsum2k dres2h = tsum2h wrk1(l1,1) = wrk1(l1,1) + dres1k * basis1 wrk1(l2,1) = wrk1(l2,1) + dres1h * basis1 wrk1(l3,1) = wrk1(l3,1) + dres1k * basis2 wrk1(l4,1) = wrk1(l4,1) + dres1h * basis2 wrk1(l5,1) = wrk1(l5,1) + dres1k * basis3 wrk1(l6,1) = wrk1(l6,1) + dres1h * basis3 wrk1(l7,1) = wrk1(l7,1) + dres1k * basis4 wrk1(l8,1) = wrk1(l8,1) + dres1h * basis4 wrk2(l1,1) = wrk2(l1,1) + dres2k * basis1 wrk2(l2,1) = wrk2(l2,1) + dres2h * basis1 wrk2(l3,1) = wrk2(l3,1) + dres2k * basis2 wrk2(l4,1) = wrk2(l4,1) + dres2h * basis2 wrk2(l5,1) = wrk2(l5,1) + dres2k * basis3 wrk2(l6,1) = wrk2(l6,1) + dres2h * basis3 wrk2(l7,1) = wrk2(l7,1) + dres2k * basis4 wrk2(l8,1) = wrk2(l8,1) + dres2h * basis4 ! ! Add the residual x quadrature weights x basis functions to r(.). ! res1 = sum1 - 1.0 res2 = sum2 - 1.0 r(l1,1) = r(l1,1) + res1 * wgtx * wgty * basis1 r(l2,1) = r(l2,1) + res2 * wgtx * wgty * basis1 r(l3,1) = r(l3,1) + res1 * wgtx * wgty * basis2 r(l4,1) = r(l4,1) + res2 * wgtx * wgty * basis2 r(l5,1) = r(l5,1) + res1 * wgtx * wgty * basis3 r(l6,1) = r(l6,1) + res2 * wgtx * wgty * basis3 r(l7,1) = r(l7,1) + res1 * wgtx * wgty * basis4 r(l8,1) = r(l8,1) + res2 * wgtx * wgty * basis4 ! ! Add derivatives of the residual x quadrature weights x basis ! functions to dr(.). ! dr(l1,1:na) = dr(l1,1:na) + wrk1(1:na,1)*wgtx*wgty*basis1 dr(l2,1:na) = dr(l2,1:na) + wrk2(1:na,1)*wgtx*wgty*basis1 dr(l3,1:na) = dr(l3,1:na) + wrk1(1:na,1)*wgtx*wgty*basis2 dr(l4,1:na) = dr(l4,1:na) + wrk2(1:na,1)*wgtx*wgty*basis2 dr(l5,1:na) = dr(l5,1:na) + wrk1(1:na,1)*wgtx*wgty*basis3 dr(l6,1:na) = dr(l6,1:na) + wrk2(1:na,1)*wgtx*wgty*basis3 dr(l7,1:na) = dr(l7,1:na) + wrk1(1:na,1)*wgtx*wgty*basis4 dr(l8,1:na) = dr(l8,1:na) + wrk2(1:na,1)*wgtx*wgty*basis4 end do end do end do end do ! ! Invert the Jacobian matrix dr ! sol = r call dgesv(na,1,dr,na,iwrk,sol,na,info) if (info/=0) write(*,*) 'WARNING: illegal value for DGESV' ! ! Do Newton update. ! do j=1,ns do i=1,nxy k = (j-1)*nxy*2+2*i-1 kbf1(i,j) = kbf0(i,j) - sol(k) hbf1(i,j) = hbf0(i,j) - sol(k+1) end do end do sum1 = dot_product(sol(:,1),sol(:,1)) sum1 = sqrt(sum1)/float(na) write(*,*) ' Residual = ',sum1 ! ! Check to see if the solution is converged. If so, stop. ! if (sum1 < crit) exit newton kbf0 = kbf1 hbf0 = hbf1 end do newton end if ! ! Simulate time series. ! if (isimul == 1) then write(*,*) write(*,*) 'Simulating time series ... ' if (idraw == 1) then ! ! Pick seed for random number generator ! useed = ustart(iseed) do j=1,nsim ! ! Randomly pick a number between 1 and ns ! rn = float(ns)*uni() do m=ns,1,-1 if (rn <= float(l)) l = m end do ind(j,1) = l ! ! Let endogenous states be steady states associated with the ! initial exogenous state. ! xt = xa(nx/2) yt = ya(ny/2) do i=1,20 x = xt y = yt ix = 1 do m = 1,nx-1 if (x > xa(m)) ix = m end do x1 = xa(ix) x2 = xa(ix+1) dx = x2-x1 iy = 1 do m = 1,ny-1 if (y > ya(m)) iy = m end do y1 = ya(iy) y2 = ya(iy+1) dy = y2-y1 nd1 = ix+nx*(iy-1) nd2 = nd1+1 nd4 = nd1+nx nd3 = nd4+1 basis1 = (x2-x)*(y2-y)/(dx*dy) basis2 = (x-x1)*(y2-y)/(dx*dy) basis3 = (x-x1)*(y-y1)/(dx*dy) basis4 = (x2-x)*(y-y1)/(dx*dy) kbf = kbf1(nd1,l)*basis1 + kbf1(nd2,l)*basis2+ & kbf1(nd3,l)*basis3 + kbf1(nd4,l)*basis4 hbf = hbf1(nd1,l)*basis1 + hbf1(nd2,l)*basis2+ & hbf1(nd3,l)*basis3 + hbf1(nd4,l)*basis4 kb = x ki = y kbi = kb-kbf yb = kbf**thet1*ki**phi1*(z1(l)*hbf)**lshare1 qxi = thet1*yb*kbi/(thet2*kbf) hbi = lshare2*qxi*hbf/(lshare1*yb) xi = kbi**thet2*ki**phi2*(z2(l)*hbi)**lshare2 q = qxi/xi h = hbf+hbi+hn(l) c = (1.0-tauh(l))*lshare1*yb*(1.0-h)/(psi*(1.0+tauc(l))*hbf) xb = yb+yn(l)-xn(l)-g(l)-c xt = ((1.0-deltb)*x+xb)/grate yt = ((1.0-delti)*y+xi)/grate end do x0(j) = xt y0(j) = yt do i=2,nt rn = uni() sum1 = 1.0 do m=ns,1,-1 if (rn <= sum1) l = m sum1 = sum1 - pi(ind(j,i-1),m) end do ind(j,i) = l end do end do endif do j=1,nsim xt = x0(j) yt = y0(j) do i=1,nt ! ! Initialize states and basis functions. ! x = xt y = yt l = ind(j,i) ix = 1 do m = 1,nx-1 if (x > xa(m)) ix = m end do x1 = xa(ix) x2 = xa(ix+1) dx = x2-x1 iy = 1 do m = 1,ny-1 if (y > ya(m)) iy = m end do y1 = ya(iy) y2 = ya(iy+1) dy = y2-y1 nd1 = ix+nx*(iy-1) nd2 = nd1+1 nd4 = nd1+nx nd3 = nd4+1 basis1 = (x2-x)*(y2-y)/(dx*dy) basis2 = (x-x1)*(y2-y)/(dx*dy) basis3 = (x-x1)*(y-y1)/(dx*dy) basis4 = (x2-x)*(y-y1)/(dx*dy) kbf = kbf1(nd1,l)*basis1 + kbf1(nd2,l)*basis2+ & kbf1(nd3,l)*basis3 + kbf1(nd4,l)*basis4 hbf = hbf1(nd1,l)*basis1 + hbf1(nd2,l)*basis2+ & hbf1(nd3,l)*basis3 + hbf1(nd4,l)*basis4 kb = x ki = y kbi = kb-kbf yb = kbf**thet1*ki**phi1*(z1(l)*hbf)**lshare1 qxi = thet1*yb*kbi/(thet2*kbf) hbi = lshare2*qxi*hbf/(lshare1*yb) xi = kbi**thet2*ki**phi2*(z2(l)*hbi)**lshare2 q = qxi/xi h = hbf+hbi+hn(l) c = (1.0-tauh(l))*lshare1*yb*(1.0-h)/(psi*(1.0+tauc(l))*hbf) xb = yb+yn(l)-xn(l)-g(l)-c xt = ((1.0-deltb)*x+xb)/grate yt = ((1.0-delti)*y+xi)/grate d = (thet1+phi1)*yb-xb-qxi-taup(l)*((thet1+phi1)*yb- & deltb*kb-tauk(l)*kb-qxi)-tauk(l)*kb+taux(l)*xb ! ! Print results for period i, simulation j. ! write(8,'(1X,I2,1X,30F11.4)') l,x,y,kbf,kbi,hbf,hbi,yb,c,h,xb, & xi,q,qxi,d,thet1*yb/kbf,phi1*yb/ki,lshare1*yb/hbf, & z1(l),z2(l),tauc(l),tauh(l),taux(l),tauk(l),taup(l), & taud(l),g(l),yn(l),xn(l),hn(l) end do write(8,*) end do endif ! ! Write out new input file with results. ! write(7,'(1X,F7.4,12X,''beta in [0,1] '')') beta write(7,'(1X,F7.4,12X,''deltb in [0,1] '')') deltb write(7,'(1X,F7.4,12X,''delti in [0,1] '')') delti write(7,'(1X,F7.4,12X,''eta >= 0 '')') eta write(7,'(1X,F7.4,12X,''gamma >= 0 '')') gamma write(7,'(1X,F7.4,12X,''thet1 in [0,1] '')') thet1 write(7,'(1X,F7.4,12X,''phi1 in [0,1] '')') phi1 write(7,'(1X,F7.4,12X,''thet2 in [0,1] '')') thet2 write(7,'(1X,F7.4,12X,''phi2 in [0,1] '')') phi2 write(7,'(1X,F7.4,12X,''psi >= 0 '')') psi write(7,'(1X,F7.4,12X,''chi in [0,1] '')') chi write(7,*) write(7,'(1X,I2,17X,''#indices in z1,z2,tauc,tauh,taux,tauk,'', & ''taup,taud,g,yn,xn,hn'')') ns write(7,*) write(7,'(30(1X,F7.4))') (z1(i), i=1,ns) write(7,'(30(1X,F7.4))') (z2(i), i=1,ns) write(7,'(30(1X,F7.4))') (tauc(i), i=1,ns) write(7,'(30(1X,F7.4))') (tauh(i), i=1,ns) write(7,'(30(1X,F7.4))') (taux(i), i=1,ns) write(7,'(30(1X,F7.4))') (tauk(i), i=1,ns) write(7,'(30(1X,F7.4))') (taup(i), i=1,ns) write(7,'(30(1X,F7.4))') (taud(i), i=1,ns) write(7,'(30(1X,F7.4))') (g(i), i=1,ns) write(7,'(30(1X,F7.4))') (yn(i), i=1,ns) write(7,'(30(1X,F7.4))') (xn(i), i=1,ns) write(7,'(30(1X,F7.4))') (hn(i), i=1,ns) write(7,*) do i=1,ns write(7,'(30(1X,F7.4))') (pi(j,i),j=1,ns) end do write(7,*) do i=1,nx-1 write(7,'(1X,F8.5)') xa(i) end do write(7,'(1X,F8.5,11X,''grid on tangible business capital'')') xa(nx) write(7,*) do i=1,ny-1 write(7,'(1X,F8.5)') ya(i) end do write(7,'(1X,F8.5,11X,''grid on intangible business capital'')') ya(ny) write(7,*) do i=1,nxy write(7,'(30(1X,F8.5))') (kbf1(i,j),j=1,ns) end do write(7,*) do i=1,nxy write(7,'(30(1X,F8.5))') (hbf1(i,j),j=1,ns) end do write(7,*) write(7,'(1X,I6,13X,''Compute equilibrium? (1=yes,0=no)'')') icomp write(7,'(1X,I6,13X,''Simulate time series? (1=yes,0=no)'')') isimul if (isimul == 1) then write(7,'(1X,I6,13X,''Number of simulations'')') nsim write(7,'(1X,I6,13X,''Number of time periods per simulation'')') nt write(7,'(1X,I6,13X,''Draw random realizations? (1=yes,0=no)'')') idraw write(7,'(1X,I6,13X,''Seed for random numbers? (integer)'')') iseed if (idraw == 0) then do i=1,nsim write(7,*) write(7,'(1X,2F8.5,70(1X,I2))') x0(i),y0(i),(ind(i,j),j=1,nt) end do endif endif ! ! Write out decision rules. ! m = 1 do j=1,ny do i=1,nx write(9,'40(1X,F8.5)') xa(i),ya(j),(kbf1(m,l),l=1,ns) m = m+1 enddo enddo write(9,*) m = 1 do j=1,ny do i=1,nx write(9,'40(1X,F8.5)') xa(i),ya(j),(hbf1(m,l),l=1,ns) m = m+1 enddo enddo end program intang