program intang ! ! DESCRIPTION: This program computes an equilibrium for the ! stochastic growth model with intangible capital and ! variation in tax rates, government consumption, and ! productivity. See great.tex and comments below for ! details. ! ! MODEL Households solve: ! ECONOMY: ! max E sum_t beta^t U(c[t],1-h[t]) N[t] ! {c,xt,xi,h} ! ! subject to ! ! c[t] + xt[t] + q[t] xi[t] ! <= rt[t] kt[t] + ri[t] ki[t] + w[t] h[t] ! - tauc[t] c[t] ! - tauh[t] w[t] h[t] ! + taub[t] (1-chi)q[t]xI[t] ! - taux[t] xt[t] ! - tauk[t] k[t] ! - taup[t]{(rt[t]-delt-tauk[t])kt[t] ! +ri[t]ki[t]-chi q[t]xi[t]} ! - tauu[t]{k[t+1]-k[t]} ! - taud[t]{(rt[t]-tauk[t]) kt[t] ! -(1+taux[t])xt[t] ! - taup[t][(rt[t]-delt-tauk[t])kt[t] ! +ri[t]ki[t]-chi q[t]xi[t]]} ! + kappa[t] ! ! kt[t+1] = [(1-delt) kt[t] + xt[t]]/(1+gn) ! ki[t+1] = [(1-deli) ki[t] + xi[t]]/(1+gn) ! ! ! Rental and Wage rates: marginal products from ! ! y[t] = F(k1t[t],ki[t],h1[t]) ! xi[t] = G(k2t[t],ki[t],h2[t]) ! ! Resource Constraint: ! ! c[t]+xt[t]+xn[t]+g[t]=y[t]+yn[t] ! ! Processes for exogenous states: ! ! i is the index for a Markov chain with transition pi ! i is in the set {1,...,ns} ! ! detrended G, taxes, technologies, and nonbusiness ! quantities are functions of i. ! ! SOLUTION Represent c and xt as weighted sums of bilinear ! basis functions: ! ! c(x,y,i) = sum_e alphac_k^i N_e(x,y) ! xt(x,y,i) = sum_e alphax_k^i N_e(x,y) ! ! where ! ! N_e(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 -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 -y_{b-1}) ! -------------* ------------- ! (x_{a+1}-x_a) (y_b-y_{b-1}) ! ! on [x_a,x_{a+1}] 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 e is at the center of the four squares in ! which N_e(x,y) is nonzero, i.e.: ! ! --*-------------*-------------* ! y_{b+1}| | | ! | | | ! | | | ! | | | ! --*-------------e-------------* ! y_b| | | ! | | | ! | | | ! | | | ! --*-------------*-------------* ! y_{b-1},x_{a-1} |x_a |x_{a+1} ! ! Plug approximate functions c(x,y,i) and xt(x,y,i) into ! the Euler equations and set them equal to 0. Assume ! that these equations can be represented as ! ! R(x,y,i;alpha) = 0. ! ! Then the approximate solutions computed here are ! the functions given above with the alpha's chosen to ! satisfy the integral equations: ! ! int R(x,y,i;alpha) N_j(x,y) dx dy = 0 ! x,y ! ! for all i,j ! ! FUNCTIONAL U(c,1-h) = log(c) + psi(1-h)^\phi/phi ! FORMS: F(k1t,ki,h1) = k1t^thet1 *ki^alph1 * h1^(1-thet1-alph1) ! G(k2t,ki,h2) = k2t^thet2 *ki^alph2 * h2^(1-thet2-alph2) ! penalty = zeta * [min(xt,0)^3 + min(xi,0)^3] ! ! TEST CASE: phi=0, delt=deli = 1, alph1=0, z2=0, all shocks off ! ! INPUT FILE: mpmodel.inp assumes the following order for inputs: ! ! alph1 : share of intangible capital in F(.) ! alph2 : share of intangible capital in G(.) ! beta : discount factor ! chi : fraction of intangible expensed (vs. sweat) ! deli : depreciation of intangible business capital ! delt : depreciation of tangible business capital ! gn : growth of population ! gz : growth of technology ! thet1 : share of tangible capital in F(.) ! thet2 : share of tangible capital in G(.) ! phi : preference parameter, curvature ! psi : preference parameter, weight on leisure ! zeta : penalty function parameter ! ns : number of exogenous states ! g : government spending (1 x ns) ! tauc : tax rate on consumption (1 x ns) ! tauh : tax rate on labor (1 x ns) ! taub : tax rate on owner's labor (1 x ns) ! taux : tax rate on on investment (1 x ns) ! tauk : tax rate on property (1 x ns) ! taup : tax rate on profits (1 x ns) ! tauu : tax rate on undistributed profits (1 x ns) ! taud : tax rate on dividends (1 x ns) ! z1 : level of technology, sector 1 (1 x ns) ! z2 : level of technology, sector 2 (1 x ns) ! hn : nonbusiness hours (1 x ns) ! xn : nonbusiness investment (1 x ns) ! yn : nonbusiness output (1 x ns) ! pi' : transpose transition matrix (ns x ns) ! (columns sum to 1) ! kt : grid on tangible capital (nkt x 1) ! ki : grid on intangible capital (nki x 1) ! init : 1=initial guesses of c,xt included, ! 0 otherwise ! c : if init=1, guess for c (nkt*nki*ns x 1) ! xt : if init=1, guess for xt (nkt*nki*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 ! kt0,ki0, : if idraw=0, initial state and indices for ! ind exogenous states (1 x 2+nt) ! ! Ellen McGrattan, 3-15-06 ! Revised, ERM, 3-17-06 implicit none !integer, parameter :: nkt=10,nki=10,na=100,ns=45,nas=9000, & integer, parameter :: nkt=10,nki=10,na=100,ns=29,nas=5800, & nq=3,msim=100,mt=2000 real, parameter :: crit= 0.000000000001 real, parameter :: crit2= 0.0000000001 real, dimension (nkt-1,nq) :: abst,wgtt real, dimension (nki-1,nq) :: absi,wgti real, dimension (nkt) :: ktn real, dimension (nki) :: kin real, dimension (ns) :: g,tauc,tauh,taub,taux,tauk,taup,tauu, & taud,z1,z2,hn,xn,yn real, dimension (ns,ns) :: pi = 0.0 real, dimension (na,ns) :: c0,c1,xt0,xt1,xi1,h11,h21,k1t1,k2t1 real, dimension (nas,1) :: r,wrk1,wrk2,rf real, dimension (nas,nas) :: dr real, dimension (msim) :: kt0,ki0 real :: alph1,alph2,beta,betg,chi,delt,deli,gn, & gz,psi,phi,thet1,thet2,zeta,delt1,deli1, & phi1,grate,lshare1,lshare2,kt1,kt2,dkt, & ki1,ki2,dki,wt,wi,wgt,c,xt,xi,k1t,k2t,kt, & k1i,k2i,ki,ct,xtt,xit,k1tt,k2tt,ktt,k1it, & k2it,kit,kt1t,kt2t,dktt,ki1t,ki2t,dkit,h, & h1,h2,h10,y,mu,dmu,ht,h1t,h2t,yt,mut, & dmut,taxc1,taxh1,taxx1,taxd1,taxc2,taxh2, & taxx2,taxp2,taxd2,chifac1,chifac2,wage, & qxi,q,rtt,rit,pent,peni,pentt,penit,rt, & ri,waget,qxit,qt,v,d,tctxt,tctktt,tctkit, & txtktt,txttkit,drttk1t,drttki,drtthi, & tctc,drtth1t,drtth1,txttktt,drttktt, & tkttxt,th1txt,txttxt,tqc,txic,tqxt,txixt, & th1tc,txttc,coef1,coef2,coef3,coefh1, & coefc,coefxt,coefkt,coefki,basis1,basis2, & basis3,basis4,basis1t,basis2t,basis3t, & basis4t,dh,dh1,dh2,dk1t,yrhs,res,dres, & dht,dh1t,dh2t,dk1tt,dresh1,dresc,dresxt, & th1c,th1xt,tyc,twc,sum1,sum2,tsum1c, & tsum2c,tsum1xt,tsum2xt,tsum1ct,tsum2ct, & tsum1xtt,tsum2xtt,tyxt,twxt,thc,thxt, & th2c,th2xt,tk1tc,tk1txt,tqxic,tqxixt, & tkic,tkixt,tqtc,tqtxt,tmuc,tkttc,tkitc, & tkitxt,rhs1,rhs2,tytc,tytct,tytxt,tytxtt, & th1tct,th1txtt,tem,thtc,thtxt,thtct, & thtxtt,th2tc,th2txt,th2tct,th2txtt, & tk1ttc,tk1ttct,tk1ttxt,tk1ttxtt,th2xtt, & tmutc,tmutxt,tmutct,tmutxtt,trttc,trttct, & trttxt,trttxtt,trhs1c,trhs1ct,trhs1xt, & trhs1xtt,trhs2c,tritc,trhs2ct,tritct, & tqtct,trhs2xt,tritxt,trhs2xtt,tritxtt, & tqtxtt,tmuxt,tk2ttc,tk2ttct,tkttct, & tk2ttxt,tk2ttxtt,tkttxtt,tqxitc,tqxitct, & tqxitxt,tqxitxtt,txitc,txitct,tkitct, & txitxt,txitxtt,tkitxtt,res1,res2,dres1c, & dres2c,dres1xt,dres2xt,uni,ustart,useed, & rn,eps integer, dimension (nas,1) :: iwrk integer, dimension (msim,mt) :: ind integer :: init,info,i,j,l,n,ne, & nd1,nd2,nd3,nd4,nd1t,nd2t,nd3t,nd4t, & s,s1,s2,s3,s4,s5,s6,s7,s8, & t,t1,t2,t3,t4,t5,t6,t7,t8, & iq,jq, & icomp,isimul,idraw,nsim,nt,iseed,it,maxit ! ! ! Parameters: ! open(unit=5, file='model3.inp') open(unit=7, file='model3.nxt') open(unit=8, file='model3.dat') open(unit=9, file='model3f.dat') open(unit=10, file='model3p.dat') read(5,*) alph1 read(5,*) alph2 read(5,*) beta read(5,*) chi read(5,*) deli read(5,*) delt read(5,*) gn read(5,*) gz read(5,*) phi read(5,*) psi read(5,*) thet1 read(5,*) thet2 read(5,*) zeta read(5,*) i if (i /= ns) then write(*,*) 'ERROR: wrong number of exogenous states. Edit mpmodels.f90' write(*,*) ' and change parameter ns.' stop endif read(5,*) g,tauc,tauh,taub,taux,tauk,taup,tauu,taud,z1,z2,hn,xn,yn read(5,*) pi read(5,*) ktn read(5,*) kin read(5,*) init if (init == 1) then do i=1,na read(5,*) (c0(i,j),j=1,ns) enddo do i=1,na read(5,*) (xt0(i,j),j=1,ns) enddo else do n=1,ns l = 1 do j=1,nki do i=1,nkt c0(l,n) = (1.0-beta*thet1)*ktn(i)**thet1*z1(n)**(1.-thet1) xt0(l,n) = beta*thet1*ktn(i)**thet1*z1(n)**(1.-thet1) l = l+1 enddo enddo enddo endif 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,*) kt0(i),ki0(i),(ind(i,j),j=1,nt) end do endif endif ! ! Intermediate parameters: ! betg = beta/(1.+gz) grate = (1.+gn)*(1.+gz) phi1 = 1.-phi h10 = 0.15 lshare1 = 1.-thet1-alph1 lshare2 = 1.-thet2-alph2 delt1 = 1.-delt deli1 = 1.-deli ! ! Quadarature abscissas and weights: ! do i=1,nkt-1 call qgausl (nq,ktn(i),ktn(i+1),abst(i,:),wgtt(i,:)) enddo do i=1,nki-1 call qgausl (nq,kin(i),kin(i+1),absi(i,:),wgti(i,:)) enddo ! ! Compute equilibrium ! if (icomp == 1) then ! ! Solve the fixed point problem: int R(x,y,z,i;phi) dx dy dz = 0 ! for phi using a Newton method. ! write(*,*) 'Computing decision functions .... ' ne = (nkt-1)*(nki-1) maxit = 10 it = 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 s=1,ns ! ! Tax coefficients ! taxc1 = 1.+tauc(s) taxx1 = 1.+taux(s) taxd1 = 1.-taud(s) taxh1 = 1.-tauh(s) chifac1 = (1.-chi)*(1.-taub(s))+chi*taxd1*(1.-taup(s)) do n=1,ne ! ! Element n is square [ktn(i),ktn(i+1)] x [kin(j),kin(j+1)] ! i = mod(n-1,nkt-1)+1 j = (n-1)/(nkt-1)+1 kt1 = ktn(i) kt2 = ktn(i+1) dkt = kt2-kt1 ki1 = kin(j) ki2 = kin(j+1) dki = ki2-ki1 nd1 = i+nkt*(j-1) nd2 = nd1+1 nd4 = nd1+nkt nd3 = nd4+1 ! ! Compute r at all quadrature points on the element. ! do jq=1,nq ki = absi(j,jq) wi = wgti(j,jq) do iq=1,nq kt = abst(i,iq) wt = wgtt(i,iq) wgt = wi*wt basis1 = (kt2-kt)*(ki2-ki)/(dkt*dki) basis2 = (kt-kt1)*(ki2-ki)/(dkt*dki) basis3 = (kt-kt1)*(ki-ki1)/(dkt*dki) basis4 = (kt2-kt)*(ki-ki1)/(dkt*dki) ! ! Compute decisions using finite element approximation. ! c = c0(nd1,s) *basis1 + c0(nd2,s) *basis2 + & c0(nd3,s) *basis3 + c0(nd4,s) *basis4 xt = xt0(nd1,s)*basis1 + xt0(nd2,s)*basis2 + & xt0(nd3,s)*basis3 + xt0(nd4,s)*basis4 !if (xt<0) write(*,'(7(1x,f7.4))') kt,ki,c,xt ! ! Back out hours of work in the sector 1 ! y = c+xt+g(s)+xn(s)-yn(s) coef1 = (taxh1*lshare1*y/(psi*taxc1*c))**(-1./phi1) coef2 = thet1*lshare2/(thet2*lshare1) coef3 = ki**alph1*z1(s)**lshare1 h1 = h10 inner_newton1: do h = 1.-coef1*h1**(1./phi1) h2 = h-h1-hn(s) k1t = coef2*h1*kt/(h2+coef2*h1) dh = -(1.-h)/(phi1*h1) dh1 = 1. dh2 = dh-dh1 dk1t = k1t*(1.-k1t/kt)*(dh1/h1-dh2/h2) yrhs = coef3*k1t**thet1*h1**lshare1 res = y-yrhs dres = -yrhs*(thet1*dk1t/k1t+lshare1/h1) if (abs(res/dres) <= crit) exit inner_newton1 h1 = h1-res/dres enddo inner_newton1 ! ! Back out other variables using static conditions. ! h = 1.-coef1*h1**(1./phi1) h2 = h-h1-hn(s) if ((h2<=0).and.(abs(phi1-1.)<1e-5)) then h2 = crit h1 = (1.-hn(s))/(1.+coef1) endif k1t = coef2*h1*kt/(h2+coef2*h1) k2t = kt-k1t wage = lshare1*y/h1 qxi = wage*h2/lshare2 xi = k2t**thet2*ki**alph2*(z2(s)*h2)**lshare2 if ((h2<=0).or.(k2t<=0)) then write(*,*) 'Found a negative factor in sector 2' write(*,'(a21,7(1x,f7.4))') 'kt,ki,c,xt,h1,h2,k2t=', & kt,ki,c,xt,h1,h2,k2t stop endif q = qxi/xi mu = 1./(taxc1*c) ktt = ((1.0-delt)*kt+xt)/grate kit = ((1.0-deli)*ki+xi)/grate pent = min(xt,0.) peni = min(xi,0.) ! ! Find tomorrow's state on the grid ! i = 1 do l = 1,nkt-1 if (ktt > ktn(l)) i = l enddo kt1t = ktn(i) kt2t = ktn(i+1) dktt = kt2t-kt1t j = 1 do l = 1,nki-1 if (kit > kin(l)) j = l enddo ki1t = kin(j) ki2t = kin(j+1) dkit = ki2t-ki1t nd1t = i+nkt*(j-1) nd2t = nd1t+1 nd4t = nd1t+nkt nd3t = nd4t+1 basis1t = (kt2t-ktt)*(ki2t-kit)/(dktt*dkit) basis2t = (ktt-kt1t)*(ki2t-kit)/(dktt*dkit) basis3t = (ktt-kt1t)*(kit-ki1t)/(dktt*dkit) basis4t = (kt2t-ktt)*(kit-ki1t)/(dktt*dkit) ! ! Compute derivatives to be used later. ! tyc = 1. tyxt = 1. coefh1 = (lshare1*h2+thet1*(1.-k1t/kt)*(h1+h2+(1.-h)/phi1))* & y/(h1*h2) coefc = 1.+thet1*(1.-k1t/kt)*(1.-h)*(1.-y/c)/(phi1*h2) coefxt = 1.+thet1*(1.-k1t/kt)*(1.-h)/(phi1*h2) th1c = coefc/coefh1 th1xt = coefxt/coefh1 twc = wage*(tyc/y-th1c/h1) twxt = wage*(tyxt/y-th1xt/h1) thc = (1.-h)/phi1*(twc/wage-1./c) thxt = (1.-h)/phi1*(twxt/wage) th2c = thc-th1c th2xt = thxt-th1xt if ((h2<=0).and.(abs(phi1-1.)<1e-5)) then write(*,*) 'in 1' th2c = 0. th2xt = 0. th1c = -h1*(coef1/c-coef1*tyc/y)/(1.+coef1) th1xt = -h1*(-coef1*tyxt/y)/(1.+coef1) endif tk1tc = k1t*(th1c/h1)-k1t*(th2c+coef2*th1c)/(h2+coef2*h1) tk1txt = k1t*(th1xt/h1)-k1t*(th2xt+coef2*th1xt)/(h2+coef2*h1) tqxic = qxi*(tyc/y+th2c/h2-th1c/h1) tqxixt = qxi*(tyxt/y+th2xt/h2-th1xt/h1) txic = xi*(-thet2*tk1tc/k2t+alph2*tkic/ki+lshare2*th2c/h2) txixt = xi*(-thet2*tk1txt/k2t+alph2*tkixt/ki+lshare2*th2xt/h2) tqc = q*(tqxic/qxi-txic/xi) tqxt = q*(tqxixt/qxi-txixt/xi) tmuc = -mu/c tkttc = 0. tkttxt = 1./grate tkitc = txic/grate tkitxt = txixt/grate ! ! Initialize the work vectors. ! wrk1 = 0. wrk2 = 0. sum1 = 0. sum2 = 0. tsum1c = 0. tsum2c = 0. tsum1xt = 0. tsum2xt = 0. ! ! Compute sums: ! ! (1) betahat sum_j pi(i,j) * ! [(1-taud')*mu' *{(1-taup')*(rt'-tauk') ! +(1-delt)*(1+taux')+delta*taup'+tauu'} ! -zeta*(1-delt)*min(xt',0)^2]/mu ! ! (2) betahat sum_j pi(i,j) * ! [mu' *{(1-taup')*(1-taud')*ri'+(1-delti)*q'* ! ((1-chi)*(1-tauh')+chi*(1-taup')*(1-taud'))} ! -zeta*(1-deli)*min(xi',0)^2]/mu ! ! do t=1,ns taxc2 = 1.+tauc(t) taxx2 = 1.+taux(t) taxd2 = 1.-taud(t) taxp2 = 1.-taup(t) taxh2 = 1.-tauh(t) chifac2 = (1.-chi)*(1.-taub(t))+chi*taxd2*taxp2 ! ! Next period consumption using finite element approximation. ! ct = c0(nd1t,t) *basis1t + c0(nd2t,t) *basis2t + & c0(nd3t,t) *basis3t + c0(nd4t,t) *basis4t xtt = xt0(nd1t,t)*basis1t + xt0(nd2t,t)*basis2t + & xt0(nd3t,t)*basis3t + xt0(nd4t,t)*basis4t ! ! Derivatives of variables with respect to coefficients ! in finite element approximation needed. ! tctktt = -(c0(nd1t,t)*basis1t+ & c0(nd4t,t)*basis4t)/(kt2t-ktt)+ & (c0(nd2t,t)*basis2t+ & c0(nd3t,t)*basis3t)/(ktt-kt1t) tctkit = -(c0(nd1t,t)*basis1t+ & c0(nd2t,t)*basis2t)/(ki2t-kit)+ & (c0(nd3t,t)*basis3t+ & c0(nd4t,t)*basis4t)/(kit-ki1t) txttktt = -(xt0(nd1t,t)*basis1t+ & xt0(nd4t,t)*basis4t)/(kt2t-ktt)+ & (xt0(nd2t,t)*basis2t+ & xt0(nd3t,t)*basis3t)/(ktt-kt1t) txttkit = -(xt0(nd1t,t)*basis1t+ & xt0(nd2t,t)*basis2t)/(ki2t-kit)+ & (xt0(nd3t,t)*basis3t+ & xt0(nd4t,t)*basis4t)/(kit-ki1t) tctc = tctktt *tkttc +tctkit *tkitc txttc = txttktt*tkttc +txttkit*tkitc tctxt = tctktt *tkttxt+tctkit *tkitxt txttxt = txttktt*tkttxt+txttkit*tkitxt ! ! Back out next period hours of work. ! yt = ct+xtt+g(t)+xn(t)-yn(t) coef1 = (taxh2*lshare1*yt/(psi*taxc2*ct))**(-1./phi1) coef2 = thet1*lshare2/(thet2*lshare1) coef3 = kit**alph1*z1(t)**lshare1 h1t = h10 inner_newton2: do ht = 1.-coef1*h1t**(1./phi1) h2t = ht-h1t-hn(t) k1tt = coef2*h1t*ktt/(h2t+coef2*h1t) dht = -(1-ht)/(phi1*h1t) dh1t = 1. dh2t = dht-dh1t dk1tt = k1tt*(1.-k1tt/ktt)*(dh1t/h1t-dh2t/h2t) yrhs = coef3*k1tt**thet1*h1t**lshare1 res = yt-yrhs dres = -yrhs*(thet1*dk1tt/k1tt+lshare1/h1t) if (abs(res/dres) <= crit) exit inner_newton2 h1t = h1t-res/dres enddo inner_newton2 ! ! Back out other variables. ! ht = 1.-coef1*h1t**(1./phi1) h2t = ht-h1t-hn(t) if ((h2t<=0).and.(abs(phi1-1.)<1e-5)) then h2t = crit h1t = (1.-hn(t))/(1.+coef1) endif k1tt = coef2*h1t*ktt/(h2t+coef2*h1t) k2tt = ktt-k1tt waget = lshare1*yt/h1t qxit = waget*h2t/lshare2 xit = k2tt**thet2*kit**alph2*(z2(t)*h2t)**lshare2 qt = qxit/xit mut = 1./(taxc2*ct) rtt = thet1*yt/k1tt rit = (alph1*yt+alph2*qxit)/kit pentt = min(xtt,0.) penit = min(xit,0.) rhs1 = taxd2*(taxp2*(rtt-tauk(t))+delt1*taxx2+delt* & taup(t)+tauu(t)) rhs2 = taxd2*taxp2*rit+deli1*qt*chifac2 ! ! Add terms to sums in Euler equations for capital. ! sum1 = sum1 + betg*pi(s,t)*(rhs1*mut-zeta*delt1* & pentt*pentt)/mu sum2 = sum2 + betg*pi(s,t)*(rhs2*mut-zeta*deli1* & penit*penit)/mu ! ! Differentiate the residual equations with respect to ! coefficients in the finite element approximation. ! tytc = tctc+txttc tytct = 1. tytxt = tctxt+txttxt tytxtt = 1. coefh1 = (lshare1*h2t+thet1*(1.-k1tt/ktt)*(h1t+h2t+ & (1.-ht)/phi1))*yt/(h1t*h2t) coefc = 1.+thet1*(1.-k1tt/ktt)*(1.-ht)*(1.-yt/ct)/(phi1*h2t) coefxt = 1.+thet1*(1.-k1tt/ktt)*(1.-ht)/(phi1*h2t) coefkt = -thet1*yt/ktt coefki = -alph1*yt/kit th1tc = (coefc*tctc+coefxt*txttc+coefkt*tkttc+ & coefki*tkitc)/coefh1 th1txt = (coefc*tctxt+coefxt*txttxt+coefkt*tkttxt+ & coefki*tkitxt)/coefh1 th1tct = coefc/coefh1 th1txtt = coefxt/coefh1 tem = -(taxh2*lshare1/(psi*taxc2))**(-1./phi1)* & (ct*h1t/yt)**(1./phi1) thtc = tem*(tctc/ct+th1tc/h1t-tytc/yt) thtxt = tem*(tctxt/ct+th1txt/h1t-tytxt/yt) thtct = tem*(1./ct+th1tct/h1t-tytct/yt) thtxtt = tem*(th1txtt/h1t-tytxtt/yt) th2tc = thtc-th1tc th2txt = thtxt-th1txt th2tct = thtct-th1tct th2txtt = thtxtt-th1txtt if ((h2<=0).and.(abs(phi1-1.)<1e-5)) then write(*,*) 'in 3' th2tc = 0. th2txt = 0. th2tct = 0. th2txtt = 0. th1tc = -h1t*(coef1*tctc/ct-coef1*tytc/yt)/(1.+coef1) th1txt = -h1t*(coef1*tctxt/ct-coef1*tytxt/yt)/(1.+coef1) th1tct = -h1t*(coef1/ct-coef1*tytct/yt)/(1.+coef1) th1txtt = -h1t*(-coef1*tytxtt/yt)/(1.+coef1) endif tk1ttc = k1tt*(tkttc/ktt+th1tc/h1t)-k1tt*(th2tc+ & coef2*th1tc)/(h2t+coef2*h1t) tk1ttct = k1tt*(th1tct/h1t)-k1tt*(th2tct+coef2*th1tct)/ & (h2t+coef2*h1t) tk1ttxt = k1tt*(tkttxt/ktt+th1txt/h1t)-k1tt*(th2txt+ & coef2*th1txt)/(h2t+coef2*h1t) tk1ttxtt = k1tt*(th1txtt/h1t)-k1tt*(th2txtt+coef2*th1txtt)/ & (h2t+coef2*h1t) tk2ttc = tkttc-tk1ttc tk2ttct = tkttct-tk1ttct tk2ttxt = tkttxt-tk1ttxt tk2ttxtt = tkttxtt-tk1ttxtt tqxitc = qxit*(tytc/yt+th2tc/h2t-th1tc/h1t) tqxitct = qxit*(tytct/yt+th2tct/h2t-th1tct/h1t) tqxitxt = qxit*(tytxt/yt+th2txt/h2t-th1txt/h1t) tqxitxtt = qxit*(tytxtt/yt+th2txtt/h2t-th1txtt/h1t) txitc = xit*(thet2*tk2ttc/k2tt+alph2*tkitc/kit+ & lshare2*th2tc/h2t) txitct = xit*(thet2*tk2ttct/k2tt+alph2*tkitct/kit+ & lshare2*th2tct/h2t) txitxt = xit*(thet2*tk2ttxt/k2tt+alph2*tkitxt/kit+ & lshare2*th2txt/h2t) txitxtt = xit*(thet2*tk2ttxtt/k2tt+alph2*tkitxtt/kit+ & lshare2*th2txtt/h2t) tqtc = qt*(tqxitc/qxit-txitc/xit) tqtct = qt*(tqxitct/qxit-txitct/xit) tqtxt = qt*(tqxitxt/qxit-txitxt/xit) tqtxtt = qt*(tqxitxtt/qxit-txitxtt/xit) tmutc = -mut*tctc/ct tmutxt = -mut*tctxt/ct tmutct = -mut/ct tmutxtt = 0. trttc = rtt*(tytc/yt-tk1ttc/k1tt) trttct = rtt*(tytct/yt-tk1ttct/k1tt) trttxt = rtt*(tytxt/yt-tk1ttxt/k1tt) trttxtt = rtt*(tytxtt/yt-tk1ttxtt/k1tt) rit = (alph1*yt+alph2*qxit)/kit tritc = alph1*yt/kit*(tytc/yt-tkitc/kit)+ & alph2*qxit/kit*(tqxitc/qxit-tkitc/kit) tritct = alph1*yt/kit*(tytct/yt-tkitct/kit)+ & alph2*qxit/kit*(tqxitct/qxit-tkitct/kit) tritxt = alph1*yt/kit*(tytxt/yt-tkitxt/kit)+ & alph2*qxit/kit*(tqxitxt/qxit-tkitxt/kit) tritxtt = alph1*yt/kit*(tytxtt/yt-tkitxtt/kit)+ & alph2*qxit/kit*(tqxitxtt/qxit-tkitxtt/kit) trhs1c = taxd2*taxp2*trttc trhs1ct = taxd2*taxp2*trttct trhs1xt = taxd2*taxp2*trttxt trhs1xtt = taxd2*taxp2*trttxtt trhs2c = taxd2*taxp2*tritc+deli1*tqtc*chifac2 trhs2ct = taxd2*taxp2*tritct+deli1*tqtct*chifac2 trhs2xt = taxd2*taxp2*tritxt+deli1*tqtxt*chifac2 trhs2xtt = taxd2*taxp2*tritxtt+deli1*tqtxtt*chifac2 tsum1c = tsum1c +betg*pi(s,t)*(tmutc*rhs1+mut*trhs1c- & mut*rhs1*tmuc/mu)/mu-betg*pi(s,t)*zeta*delt1* & (2.*pentt*txttc-pentt*pentt*tmuc/mu)/mu tsum1xt = tsum1xt+betg*pi(s,t)*(tmutxt*rhs1+mut*trhs1xt)/ & mu-betg*pi(s,t)*zeta*delt1*(2.*pentt*txttxt- & pentt*pentt*tmuxt/mu)/mu tsum2c = tsum2c +betg*pi(s,t)*(tmutc*rhs2+mut*trhs2c- & rhs2*mut*tmuc/mu)/mu-betg*pi(s,t)*zeta*deli1* & (2.*penit*txttc-penit*penit*tmuc/mu)/mu tsum2xt = tsum2xt+betg*pi(s,t)*(tmutxt*rhs2+mut*trhs2xt)/ & mu-betg*pi(s,t)*zeta*deli1*(2.*penit*txitxt- & penit*penit*tmuxt/mu)/mu tsum1ct = betg*pi(s,t)*(tmutct*rhs1+mut*trhs1ct)/mu tsum2ct = betg*pi(s,t)*(tmutct*rhs2+mut*trhs2ct)/mu tsum1xtt = betg*pi(s,t)*(mut*trhs1xtt-2.*zeta*delt1*pentt)/mu tsum2xtt = betg*pi(s,t)*(mut*trhs2xtt-2.*zeta*deli1*penit* & txitxtt)/mu ! ! Hold results in work vectors for later. ! t1 = 2*nd1t-1+(t-1)*2*na t2 = t1+1 t3 = t1+2 t4 = t3+1 t5 = t3+2*nkt t6 = t5+1 t7 = t1+2*nkt t8 = t7+1 wrk1(t1,1)= wrk1(t1,1) + tsum1ct * basis1t wrk1(t2,1)= wrk1(t2,1) + tsum1xtt * basis1t wrk1(t3,1)= wrk1(t3,1) + tsum1ct * basis2t wrk1(t4,1)= wrk1(t4,1) + tsum1xtt * basis2t wrk1(t5,1)= wrk1(t5,1) + tsum1ct * basis3t wrk1(t6,1)= wrk1(t6,1) + tsum1xtt * basis3t wrk1(t7,1)= wrk1(t7,1) + tsum1ct * basis4t wrk1(t8,1)= wrk1(t8,1) + tsum1xtt * basis4t wrk2(t1,1)= wrk2(t1,1) + tsum2ct * basis1t wrk2(t2,1)= wrk2(t2,1) + tsum2xtt * basis1t wrk2(t3,1)= wrk2(t3,1) + tsum2ct * basis2t wrk2(t4,1)= wrk2(t4,1) + tsum2xtt * basis2t wrk2(t5,1)= wrk2(t5,1) + tsum2ct * basis3t wrk2(t6,1)= wrk2(t6,1) + tsum2xtt * basis3t wrk2(t7,1)= wrk2(t7,1) + tsum2ct * basis4t wrk2(t8,1)= wrk2(t8,1) + tsum2xtt * basis4t enddo ! ! The vectors wrk1 and wrk2 contain derivatives of the residuals. ! s1 = 2*nd1-1+(s-1)*2*na s2 = s1+1 s3 = s1+2 s4 = s3+1 s5 = s3+2*nkt s6 = s5+1 s7 = s1+2*nkt s8 = s7+1 dres1c = tsum1c -zeta*pent*pent*tmuc/(mu*mu) dres1xt = tsum1xt + 2.0*zeta*pent/mu dres2c = tsum2c -chifac1*tqc + 2.0*zeta*peni*txic/mu - & zeta*peni*peni*tmuc/(mu*mu) dres2xt = tsum2xt -chifac1*tqxt + 2.0*zeta*peni*txixt/mu wrk1(s1,1) = wrk1(s1,1) + dres1c * basis1 wrk1(s2,1) = wrk1(s2,1) + dres1xt * basis1 wrk1(s3,1) = wrk1(s3,1) + dres1c * basis2 wrk1(s4,1) = wrk1(s4,1) + dres1xt * basis2 wrk1(s5,1) = wrk1(s5,1) + dres1c * basis3 wrk1(s6,1) = wrk1(s6,1) + dres1xt * basis3 wrk1(s7,1) = wrk1(s7,1) + dres1c * basis4 wrk1(s8,1) = wrk1(s8,1) + dres1xt * basis4 wrk2(s1,1) = wrk2(s1,1) + dres2c * basis1 wrk2(s2,1) = wrk2(s2,1) + dres2xt * basis1 wrk2(s3,1) = wrk2(s3,1) + dres2c * basis2 wrk2(s4,1) = wrk2(s4,1) + dres2xt * basis2 wrk2(s5,1) = wrk2(s5,1) + dres2c * basis3 wrk2(s6,1) = wrk2(s6,1) + dres2xt * basis3 wrk2(s7,1) = wrk2(s7,1) + dres2c * basis4 wrk2(s8,1) = wrk2(s8,1) + dres2xt * basis4 ! ! Add the residuals x quadrature weights x basis functions to r(.). ! res1 = sum1 - taxd1*(taxx1+tauu(s))+zeta*pent*pent/mu res2 = sum2 - q*chifac1 + zeta*peni*peni/mu r(s1,1) = r(s1,1) + res1 * wgt * basis1 r(s2,1) = r(s2,1) + res2 * wgt * basis1 r(s3,1) = r(s3,1) + res1 * wgt * basis2 r(s4,1) = r(s4,1) + res2 * wgt * basis2 r(s5,1) = r(s5,1) + res1 * wgt * basis3 r(s6,1) = r(s6,1) + res2 * wgt * basis3 r(s7,1) = r(s7,1) + res1 * wgt * basis4 r(s8,1) = r(s8,1) + res2 * wgt * basis4 ! ! Add derivatives of the residuals x quadrature weights x basis ! functions to dr(.). ! dr(s1,1:nas) = dr(s1,1:nas) + wrk1(1:nas,1) * wgt * basis1 dr(s2,1:nas) = dr(s2,1:nas) + wrk2(1:nas,1) * wgt * basis1 dr(s3,1:nas) = dr(s3,1:nas) + wrk1(1:nas,1) * wgt * basis2 dr(s4,1:nas) = dr(s4,1:nas) + wrk2(1:nas,1) * wgt * basis2 dr(s5,1:nas) = dr(s5,1:nas) + wrk1(1:nas,1) * wgt * basis3 dr(s6,1:nas) = dr(s6,1:nas) + wrk2(1:nas,1) * wgt * basis3 dr(s7,1:nas) = dr(s7,1:nas) + wrk1(1:nas,1) * wgt * basis4 dr(s8,1:nas) = dr(s8,1:nas) + wrk2(1:nas,1) * wgt * basis4 enddo enddo enddo enddo ! ! Solve dr\r and put solution in r. Issue warnings if necessary. ! call dgesv(nas,1,dr,nas,iwrk,r,nas,info) if (info < 0) write(*,*) 'WARNING: illegal value for', & ' argument ',-info,' in intang.f90' if (info > 0) write(*,*) 'WARNING: singular matrix', & ' when solving dr*x=r in intang.f90' ! ! Do Newton update. ! do s=1,ns do i=1,na j = (s-1)*na*2+2*i-1 c1(i,s) = c0(i,s) - r(j,1) xt1(i,s) = xt0(i,s) - r(j+1,1) enddo enddo sum1 = dot_product(r(:,1),r(:,1)) sum1 = sqrt(sum1)/float(nas) write(*,*) ' Residual = ',sum1 ! ! Check to see if the solution is converged. If so, stop. ! if ((sum1 < crit2).or.(it==maxit)) exit newton c0 = c1 xt0 = xt1 it = it+1 enddo newton endif ! ! Simulate time series. ! if (isimul == 1) then write(*,*) write(*,*) 'Simulating time series for intangsim.m ... ' 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 l=ns,1,-1 if (rn <= float(l)) s = l end do ind(j,1) = s taxc1 = 1.+tauc(s) taxh1 = 1.-tauc(s) ! ! Let capital stocks be steady states associated with the ! initial exogenous state. ! ktt = ktn(nkt/2) kit = kin(nki/2) do i=1,20 kt = ktt ki = kit iq = 1 do l = 1,nkt-1 if (kt > ktn(l)) iq = l end do kt1 = ktn(iq) kt2 = ktn(iq+1) jq = 1 do l = 1,nki-1 if (ki > kin(l)) jq = l end do ki1 = kin(jq) ki2 = kin(jq+1) dkt = kt2-kt1 dki = ki2-ki1 basis1 = (kt2-kt)*(ki2-ki)/(dkt*dki) basis2 = (kt-kt1)*(ki2-ki)/(dkt*dki) basis3 = (kt-kt1)*(ki-ki1)/(dkt*dki) basis4 = (kt2-kt)*(ki-ki1)/(dkt*dki) nd1 = iq+nkt*(jq-1) nd2 = nd1+1 nd4 = nd1+nkt nd3 = nd4+1 ! ! Compute consumption using finite element approximation. ! c = c1(nd1,s) *basis1 + c1(nd2,s) *basis2 + & c1(nd3,s) *basis3 + c1(nd4,s) *basis4 xt = xt1(nd1,s)*basis1 + xt1(nd2,s)*basis2 + & xt1(nd3,s)*basis3 + xt1(nd4,s)*basis4 ! ! Back out hours of work in the sector 1 ! y = c+xt+g(s)+xn(s)-yn(s) coef1 = (taxh1*lshare1*y/(psi*taxc1*c))**(-1./phi1) coef2 = thet1*lshare2/(thet2*lshare1) coef3 = ki**alph1*z1(s)**lshare1 h1 = h10 inner_newton3: do h = 1.-coef1*h1**(1./phi1) h2 = h-h1-hn(s) k1t = coef2*h1*kt/(h2+coef2*h1) dh = -(1-h)/(phi1*h1) dh1 = 1. dh2 = dh-dh1 dk1t = k1t*(1.-k1t/kt)*(dh1/h1-dh2/h2) yrhs = coef3*k1t**thet1*h1**lshare1 res = y-yrhs dres = -yrhs*(thet1*dk1t/k1t+lshare1/h1) if (abs(res/dres) <= crit) exit inner_newton3 h1 = h1-res/dres enddo inner_newton3 ! ! Back out other variables using static conditions. ! h = 1.-coef1*h1**(1./phi1) h2 = h-h1-hn(s) if ((h2<=0).and.(abs(phi1-1.)<1e-5)) then h2 = crit h1 = (1.-hn(s))/(1.+coef1) endif k1t = coef2*h1*kt/(h2+coef2*h1) k2t = kt-k1t wage = lshare1*y/h1 qxi = wage*h2/lshare2 xi = k2t**thet2*ki**alph2*(z2(s)*h2)**lshare2 q = qxi/xi ! ! Update capital stocks. ! ktt = ((1.0-delt)*kt+xt)/grate kit = ((1.0-deli)*ki+xi)/grate end do kt0(j) = ktt ki0(j) = kit do i=2,nt rn = uni() sum1 = 1. do l=ns,1,-1 if (rn <= sum1) s = l sum1 = sum1 - pi(ind(j,i-1),l) end do ind(j,i) = s end do end do endif do j=1,nsim ktt = kt0(j) kit = ki0(j) do i=1,nt ! ! Initialize states and basis functions. ! kt = ktt ki = kit s = ind(j,i) taxc1 = 1.+tauc(s) taxh1 = 1.-tauh(s) iq = 1 do l = 1,nkt-1 if (kt > ktn(l)) iq = l end do kt1 = ktn(iq) kt2 = ktn(iq+1) jq = 1 do l = 1,nki-1 if (ki > kin(l)) jq = l end do ki1 = kin(jq) ki2 = kin(jq+1) dkt = kt2-kt1 dki = ki2-ki1 basis1 = (kt2-kt)*(ki2-ki)/(dkt*dki) basis2 = (kt-kt1)*(ki2-ki)/(dkt*dki) basis3 = (kt-kt1)*(ki-ki1)/(dkt*dki) basis4 = (kt2-kt)*(ki-ki1)/(dkt*dki) nd1 = iq+nkt*(jq-1) nd2 = nd1+1 nd4 = nd1+nkt nd3 = nd4+1 ! ! Compute consumption using finite element approximation. ! c = c1(nd1,s) *basis1 + c1(nd2,s) *basis2 + & c1(nd3,s) *basis3 + c1(nd4,s) *basis4 xt = xt1(nd1,s)*basis1 + xt1(nd2,s)*basis2 + & xt1(nd3,s)*basis3 + xt1(nd4,s)*basis4 ! ! Back out hours of work in the sector 1 ! y = c+xt+g(s)+xn(s)-yn(s) coef1 = (taxh1*lshare1*y/(psi*taxc1*c))**(-1./phi1) coef2 = thet1*lshare2/(thet2*lshare1) coef3 = ki**alph1*z1(s)**lshare1 h1 = h10 inner_newton4: do h = 1.-coef1*h1**(1./phi1) h2 = h-h1-hn(s) k1t = coef2*h1*kt/(h2+coef2*h1) dh = -(1-h)/(phi1*h1) dh1 = 1. dh2 = dh-dh1 dk1t = k1t*(1.-k1t/kt)*(dh1/h1-dh2/h2) yrhs = coef3*k1t**thet1*h1**lshare1 res = y-yrhs dres = -yrhs*(thet1*dk1t/k1t+lshare1/h1) if (abs(res/dres) <= crit) exit inner_newton4 h1 = h1-res/dres enddo inner_newton4 ! ! Back out other variables using static conditions. ! h = 1.-coef1*h1**(1./phi1) h2 = h-h1-hn(s) if ((h2<=0).and.(abs(phi1-1.)<1e-5)) then h2 = crit h1 = (1.-hn(s))/(1.+coef1) !write(*,*) 'Bound hit at period ',i endif k1t = coef2*h1*kt/(h2+coef2*h1) k2t = kt-k1t wage = lshare1*y/h1 qxi = wage*h2/lshare2 rt = thet1*y/k1t ri = (alph1*y+alph2*qxi)/ki xi = k2t**thet2*ki**alph2*(z2(s)*h2)**lshare2 q = qxi/xi ! ! Update capital stocks. ! ktt = ((1.-delt)*kt+xt)/grate kit = ((1.-deli)*ki+xi)/grate v = (1.-taud(s))*(1.+taux(s))*(1.+tauu(s))*ktt+((1.-chi)* & (1.-taub(s))+chi*(1.-taud(s))*(1.-taup(s)))*q*kit d = rt*kt-xt+ri*ki-chi*qxi-taup(s)*((rt-delt-tauk(s))*kt+ & ri*ki-chi*qxi)-taux(s)*xt-tauk(s)*kt-tauu(s)*(xt-delt*kt) ! ! Print results for period i, simulation j. ! write(8,'(1X,I2,1X,32F12.6)') s,kt,ki,ktt,kit,c,h1,h2,xt,xi,q,y, & rt,ri,wage,g(s),tauc(s),tauh(s),taub(s),taux(s),tauk(s),taup(s), & tauu(s),taud(s),z1(s),z2(s),hn(s),xn(s),yn(s),v,d,k1t,k2t end do write(8,*) end do endif ! ! Write out new input file with results. ! write(7,'(1X,F7.4,6X,''alph1 in [0,1] '')') alph1 write(7,'(1X,F7.4,6X,''alph2 in [0,1] '')') alph2 write(7,'(1X,F7.4,6X,''beta in [0,1] '')') beta write(7,'(1X,F7.4,6X,''chi in [0,1] '')') chi write(7,'(1X,F7.4,6X,''deli in [0,1] '')') deli write(7,'(1X,F7.4,6X,''delt in [0,1] '')') delt write(7,'(1X,F7.4,6X,''gn '')') gn write(7,'(1X,F7.4,6X,''gz '')') gz write(7,'(1X,F7.4,6X,''phi <= 0 '')') phi write(7,'(1X,F7.4,6X,''psi <= 0 '')') psi write(7,'(1X,F7.4,6X,''thet1 in [0,1] '')') thet1 write(7,'(1X,F7.4,6X,''thet2 in [0,1] '')') thet2 write(7,'(1X,F7.1,6X,''zeta >= 0 '')') zeta write(10,*) alph1 write(10,*) alph2 write(10,*) beta write(10,*) chi write(10,*) deli write(10,*) delt write(10,*) gn write(10,*) gz write(10,*) phi write(10,*) psi write(10,*) thet1 write(10,*) thet2 write(7,*) write(7, & '(1X,I2,13X,''# in g,tauc,tauh,taub,taux,tauk,taup,tauu,taud,z1,z2,'', & ''hn,xn,yn,pi below'')') ns write(7,*) write(7,'(50(1X,F7.4))') (g(i), i=1,ns) write(7,'(50(1X,F7.4))') (tauc(i), i=1,ns) write(7,'(50(1X,F7.4))') (tauh(i), i=1,ns) write(7,'(50(1X,F7.4))') (taub(i), i=1,ns) write(7,'(50(1X,F7.4))') (taux(i), i=1,ns) write(7,'(50(1X,F7.4))') (tauk(i), i=1,ns) write(7,'(50(1X,F7.4))') (taup(i), i=1,ns) write(7,'(50(1X,F7.4))') (tauu(i), i=1,ns) write(7,'(50(1X,F7.4))') (taud(i), i=1,ns) write(7,'(50(1X,F7.4))') (z1(i), i=1,ns) write(7,'(50(1X,F7.4))') (z2(i), i=1,ns) write(7,'(50(1X,F7.4))') (hn(i), i=1,ns) write(7,'(50(1X,F7.4))') (xn(i), i=1,ns) write(7,'(50(1X,F7.4))') (yn(i), i=1,ns) write(7,*) do i=1,ns write(7,'(50(1X,F7.4))') (pi(j,i),j=1,ns) enddo write(7,*) do i=1,nkt-1 write(7,'(1X,F8.5)') ktn(i) enddo write(7,'(1X,F8.5,11X,''grid on tangible capital'')') ktn(nkt) do i=1,nki-1 write(7,'(1X,F8.5)') kin(i) enddo write(7,'(1X,F8.5,11X,''grid on intangible capital'')') kin(nki) write(7,*) write(7,'(2X,I1,17X,''Input initial consumption? (1=yes,0=no)'')') init write(7,*) do i=1,na write(7,'(50(1X,F8.5))') (c1(i,j),j=1,ns) enddo do i=1,na write(7,'(50(1X,F8.5))') (xt1(i,j),j=1,ns) enddo 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,'(2(1X,F8.5),70(1X,I2))') kt0(i),ki0(i),(ind(i,j),j=1,nt) enddo endif endif ! ! Write out decision rules ! do s=1,ns taxc1 = 1.+tauc(s) taxh1 = 1.-tauh(s) taxx1 = 1.+taux(s) taxd1 = 1.-taud(s) chifac1 = (1.-chi)*(1.-taub(s))+chi*taxd1*(1.-taup(s)) do j=1,nki ki = kin(j) do i=1,nkt kt = ktn(i) c = c1((j-1)*nkt+i,s) xt = xt1((j-1)*nkt+i,s) y = c+xt+g(s)+xn(s)-yn(s) coef1 = (taxh1*lshare1*y/(psi*taxc1*c))**(-1./phi1) coef2 = thet1*lshare2/(thet2*lshare1) coef3 = ki**alph1*z1(s)**lshare1 h1 = h10 inner_newton5: do h = 1.-coef1*h1**(1./phi1) h2 = h-h1-hn(s) k1t = coef2*h1*kt/(h2+coef2*h1) dh = -(1.-h)/(phi1*h1) dh1 = 1. dh2 = dh-dh1 dk1t = k1t*(1.-k1t/kt)*(dh1/h1-dh2/h2) yrhs = coef3*k1t**thet1*h1**lshare1 res = y-yrhs dres = -yrhs*(thet1*dk1t/k1t+lshare1/h1) if (abs(res/dres) <= crit) exit inner_newton5 h1 = h1-res/dres enddo inner_newton5 h = 1.-coef1*h1**(1./phi1) h2 = h-h1-hn(s) if ((h2<=0).and.(abs(phi1-1.)<1e-5)) then h2 = crit h1 = (1.-hn(s))/(1.+coef1) endif k1t = coef2*h1*kt/(h2+coef2*h1) k2t = kt-k1t wage = lshare1*y/h1 qxi = wage*h2/lshare2 xi = k2t**thet2*ki**alph2*(z2(s)*h2)**lshare2 xi1((j-1)*nkt+i,s) = xi h11((j-1)*nkt+i,s) = h1 h21((j-1)*nkt+i,s) = h2 k1t1((j-1)*nkt+i,s) = k1t k2t1((j-1)*nkt+i,s) = k2t enddo enddo enddo do i=1,na write(9,'(30(1X,F8.5))') (c1(i,j),j=1,ns) enddo do i=1,na write(9,'(30(1X,F8.5))') (xt1(i,j),j=1,ns) enddo do i=1,na write(9,'(30(1X,F8.5))') (xi1(i,j),j=1,ns) enddo do i=1,na write(9,'(30(1X,F8.5))') (h11(i,j),j=1,ns) enddo do i=1,na write(9,'(30(1X,F8.5))') (h21(i,j),j=1,ns) enddo do i=1,na write(9,'(30(1X,F8.5))') (k1t1(i,j),j=1,ns) enddo do i=1,na write(9,'(30(1X,F8.5))') (k2t1(i,j),j=1,ns) enddo end program intang