### VERSION 4.8.2008 ### # Housekeeping reset; option show_stats 1; option solver "knitroampl"; model kmax.mod; data st.dat; solve; ## OUTPUT ## printf "\n\n******* Output for kmax ********"; printf "\n K at good steady state: %f\n\n", kmax; model sttom.mod; ### Generates output matrix of variables evaluated at the Chebyshev nodes ### set Kspace:= 1..nn; param ntil {Kspace}; param cytil{Kspace}; param cogtil{Kspace}; param cobtil{Kspace}; param Xtil{Kspace}; param Itil{Kspace}; #printf "param MATRIXtarik : 1 2 3 4 5 6 7 8 9:=\n" >OUTtarik.out; for {kindex in 1..nn} { let kinit := K[kindex]; let X:=0.2; let I:=0.4; let Rb:=0.1; let Rg:=0.3; let cy:=0.2; let cog:=0.2; let cob:=0.2; let nnew:=0.5; repeat while (abs(nnew-n) > tolerance) { if (nnew < 1) then { let n:=nnew; solve; let nnew:=(I*(1-gamma)/D)+gamma; } if (nnew >=1) then { let n:=1; solve; let nnew:=I*(1-gamma)/D+gamma; if nnew > 1 then let nnew:=1; } } if n==1 then let Rb:=0; if n==1 then let cob:=0; if n==1 then let Rg:=alpha*(Q*I)^(alpha-1); if n==1 then let cy:=(1-alpha)*(kinit^alpha)-I; if n==1 then let cog:=Rg*Q*I; printf "%d %f %f %f %f %f %f %f %f %f\n", kindex, K[kindex], n, X, I, Rg, Rb, cog, cob, cy >CWeights.out; let ntil[kindex]:=n; let cytil[kindex]:=cy; let cobtil[kindex]:=cob; let cogtil[kindex]:=cog; let Xtil[kindex]:=X; let Itil[kindex]:=I; } #### Generate the Chebyshev weights param Corder := nn; #actually of order Corder -1, just need # vectors where the weights will go: param n_CWeights {k in 1..Corder}; param cy_CWeights {k in 1..Corder}; param cob_CWeights {k in 1..Corder}; param cog_CWeights {k in 1..Corder}; param X_CWeights {k in 1..Corder}; param I_CWeights {k in 1..Corder}; param Tjm{Kspace}; let n_CWeights[1]:=(1/nn)*sum{k in Kspace} ntil[k]; for {mm in 2..Corder} { for {j in Kspace}{ let Tjm[j] := cos((mm-1)*acos(Chebnodes[j])); #the Chebyshev polynomial } let n_CWeights[mm]:=2/nn*sum{x in Kspace} (Tjm[x]*ntil[x]); } let cy_CWeights[1]:=(1/nn)*sum{k in Kspace} cytil[k]; for {mm in 2..Corder} { for {j in Kspace}{ let Tjm[j] := cos((mm-1)*acos(Chebnodes[j])); #the Chebyshev polynomial } let cy_CWeights[mm]:=2/nn*sum{x in Kspace} (Tjm[x]*cytil[x]); } let cob_CWeights[1]:=(1/nn)*sum{k in Kspace} cobtil[k]; for {mm in 2..Corder} { for {j in Kspace}{ let Tjm[j] := cos((mm-1)*acos(Chebnodes[j])); #the Chebyshev polynomial } let cob_CWeights[mm]:=2/nn*sum{x in Kspace} (Tjm[x]*cobtil[x]); } let cog_CWeights[1]:=(1/nn)*sum{k in Kspace} cogtil[k]; for {mm in 2..Corder} { for {j in Kspace}{ let Tjm[j] := cos((mm-1)*acos(Chebnodes[j])); #the Chebyshev polynomial } let cog_CWeights[mm]:=2/nn*sum{x in Kspace} (Tjm[x]*cogtil[x]); } let X_CWeights[1]:=(1/nn)*sum{k in Kspace} Xtil[k]; for {mm in 2..Corder} { for {j in Kspace}{ let Tjm[j] := cos((mm-1)*acos(Chebnodes[j])); #the Chebyshev polynomial } let X_CWeights[mm]:=2/nn*sum{x in Kspace} (Tjm[x]*Xtil[x]); } let I_CWeights[1]:=(1/nn)*sum{k in Kspace} Itil[k]; for {mm in 2..Corder} { for {j in Kspace}{ let Tjm[j] := cos((mm-1)*acos(Chebnodes[j])); #the Chebyshev polynomial } let I_CWeights[mm]:=2/nn*sum{x in Kspace} (Tjm[x]*Itil[x]); } # Display weights of the Chebyshev polynomials printf "param : n_CWeights cy_CWeights cog_CWeights cob_CWeights X_CWeights I_CWeights :=\n" > "CWeights.out"; for {line in 1..nn}{ printf "%d %f %f %f %f %f %f \n", line, n_CWeights[line], cy_CWeights[line],cog_CWeights[line], cob_CWeights[line], X_CWeights[line], I_CWeights[line] > "CWeights.out"; } printf ";" > "CWeights.out"; ##################### # SIMULATIONS # ##################### ######################### # Variables definition # ######################### var k; param periods=40; #Number of periods per simulation param nval; param outtom; param meangrowth; param vargrowth; param meanoutput; param varoutput; param prob; param captom; param ytrajectory{1..periods}; param ygrowth{1..(periods-1)}; #additional variables to compute growth rates var pre ; fix pre := k; #intialize var post; # Define lower and upper bounds to do simulation param aux1; param aux2; let aux1 := 0.4*kmax; let aux2 := 0.6* kmax; # Number of simulations param numsim = 20; param varout{1..numsim}; param meanout{1..numsim}; param kvalues{1..numsim}:= Uniform(aux1,aux2); display kvalues; #let k:=kvalues[1]; #display k; # Call model and data data CWeights.out; for {v in 1..numsim} { display v; let k := kvalues[v]; let pre := kvalues[v]; display k; for {ww in 1..periods}{ # display ww; #calculate n for a given k let nval:= sum {i in 1..nn} (n_CWeights[i]*cos((i-1)*acos(2*(k-klow)/(khigh-klow)-1))); let prob := Uniform(0,1); display prob; #param prob{nval} = Uniform(0,1); # compute capital tomorrow (note that interest is multiplied) if prob<=nval then let captom := ( alpha*( sum {i in 1..nn} (q*X_CWeights[i]*cos((i-1)*acos(2*(k-klow)/(khigh-klow)-1) ) +Q*I_CWeights[i]*cos( (i-1)*acos(2*(k-klow)/(khigh-klow)-1))))^(alpha) ); else let captom := ( alpha*( sum {i in 1..nn} (q*X_CWeights[i]*cos((i-1)*acos(2*(k-klow)/(khigh-klow)-1))))^(alpha)); # Redefine capital display k; display captom; printf "\n assigning captom to k and post\n"; unfix k; fix k := captom; display k; fix post := captom; display post; # compute output tomorrow (sum young tomorrow to old tomorrow) if prob<=nval then let outtom := (sum {i in 1..nn} (cy_CWeights[i]+cog_CWeights[i])*cos((i-1)*acos((2*(k-klow)/(kmax-klow)-1))) ); else let outtom := (sum {i in 1..nn} (cy_CWeights[i]+cob_CWeights[i])*cos((i-1)*acos((2*(k-klow)/(kmax-klow)-1))) ); # store values in arrays let ytrajectory[ww]:=outtom; printf "%d %f \n", ww, ytrajectory[ww]>ytrajectory.out; #let ygrwth[ww]:=100*(post-pre)/(pre)}; if wwygrowth.out; unfix pre; fix pre := post; #store for next period unfix post; display ww; }; # Compute mean output growth let meangrowth := (sum {i in 1..(nn-1)} ygrowth[i])/(nn-1); #Compute variance of growth let vargrowth := (sum {i in 1..(nn-1)} (ygrowth[i]-meangrowth)^2)/(nn-2); # Compute variance of output let meanoutput := (sum {i in 1..nn} ytrajectory[i] )/(nn); let varoutput := (sum {i in 1..(nn-1)} (ygrowth[i]-meanoutput)^2)/(nn-1); display ygrowth; display ytrajectory; display varoutput,vargrowth,meangrowth; display varoutput; let varout[v] := varoutput; let meanout[v] := meanoutput; }; ######################## # Summary statistics # ######################## param MEAN = (sum {i in 1..numsim} meanout[i])/(numsim); param meanvar= (sum {i in 1..(numsim-1)} varout[i])/(numsim-1); param VAR = (sum {i in 1..numsim} (varout[i]-meanvar)^2)/(numsim-2); display MEAN, VAR; printf "%f %f\n", ytrajectory>ytrajectory[psi].out; printf "%f %f\n", ygrowth>ygrowth[psi].out;