I recently try to modify sms_5_2.rpf in "RATS Handbook for Switching Models and Structural Breaks" for my own use. My model has three variables. The original code only generates one S.D. shock for girf, so I first modified the code in the following way to generate one unit shock.
********************************************
do draw=1,ndraws
gset bishocks = %ranmvnormal(fsigma)
set fshocks(1) = bishocks(t)(1)
set fshocks(2) = bishocks(t)(2)
set fshocks(3) = bishocks(t)(3)
forecast(paths,model=tvecm,results=basesims)
# fshocks
compute ishock=%ran(fsigma(1,1))
compute fshocks(1)(baseentry)=bishocks(baseentry)(1)+ishock
compute fshocks(2)(baseentry)=bishocks(baseentry)(2)
compute fshocks(3)(baseentry)=bishocks(baseentry)(3)
forecast(paths,model=tvecm,results=sims)
# fshocks
do i=1,6
set girf1(i) = girf1(i)+(sims(i)-basesims(i))/ishock
*save my each draw of girf to a matrix called sd_er
set sd_er(draw,i) = (sims(i)-basesims(i))/ishock
end do i
end do draw
do i=1,6
set girf1(i) = girf1(i)/ndraws
end do i
*************************************************************************************
Then I generate my error bands using the following code. First, I calculate standard deviation for 5000 draws of girf in each period. Second, I compute upper limit by adding one S.D. to average girf and compute lower limit by subtracting one S.D. to average girf. However, it is identical to my average girfs. That is my generated bands are the same as my girfs. I doubt there is something wrong with the matrix that saves all of my draws, but I don't know how to make a correction.
***********************************************************
* 84% upper limit, 16% lower limit
*
do i=1,6
do t=baseentry,baseentry+nsteps-1
do draw=1,ndraws
compute girf1_up(i)(t)=girf1_up(i)(t)+(sd_er(draw,i)(t)-girf1(i)(t))^2
end do draw
compute girf1_low(i)(t)=girf1(i)(t)-sqrt(girf1_up(i)(t)/(ndraws-1))
compute girf1_up(i)(t)=girf1(i)(t)+sqrt(girf1_up(i)(t)/(ndraws-1))
end do t
end do i
*************************************************************************
For your convenience, I attached my full code and the dataset. Any help will be highly appreciated. Thanks a million in advance!
Best Regards,
Yu You
Code: Select all
OPEN DATA main_ca_final_rats_totreer.dta
CALENDAR(PANELOBS=15,A) 1995:1
DATA(FORMAT=DTA) 1//1995:01 58//2009:01 code_ifs year $
ka kai kao eq eqi eqo bo boi boo mm mmi mmo ci cii cio fc fci fco di dii_ldi dii dio ldi pop xrat ppp $
cgdp cc cg price_lv openc relative_inc rgdpl kc kg ki openk ers mi ka_open base cu peg fin_dvp inflation $
trade_open_wdi rgdp_wdi nation lngr pi eqci kac bofc dbacba llgdp cbagdp dbagdp ofagdp pcrdbgdp pcrdbofgdp $
bdgdp fdgdp bcbd ll_usd overhead netintmargin concentration roa roe costinc zscore inslife insnonlife $
stmktcap stvaltraded stturnover listco_pc prbond pubond intldebt intldebtnet nrbloan offdep remit eqci_peg $
bofc_peg di_peg eqcii eqcio bofci bofco eqcii_peg eqcio_peg bofci_peg bofco_peg dii_peg dio_peg ka_peg $
kai_peg kao_peg eqpeg bopeg mmpeg cipeg fcpeg dipeg eqci_ers bofc_ers di_ers eqcii_ers eqcio_ers bofci_ers $
bofco_ers dii_ers dio_ers ka_ers kai_ers kao_ers growth lnrgdp eq_peg bo_peg life_exp lit_adu gov_deb $
pri_sch_gro sec_sch_gro kaci kaco lnpop outflow outflow_di outflow_port outflow_eq outflow_debt inflow $
inflow_di inflow_port inflow_eq inflow_debt inflowg inflowg_di inflowg_port inflowg_eq inflowg_debt $
outflowg outflowg_di outflowg_port outflowg_eq outflowg_debt flowg flowg_di flowg_port flowg_eq flowg_debt $
ka_flowg kai_inflowg kao_outflowg eq_flowg eqi_inflowg eqo_outflowg debt_flowg debti_inflowg debto_outflowg $
di_flowg dii_inflowg dio_outflowg debti debto debt gov_eff gov_eff_rk reg_qua reg_qua_rk rul_law rul_law_rk $
con_cor con_cor_rk voi_acc voi_acc_rk pol_sta pol_sta_rk gov_eff_ka rul_law_ka con_cor_ka reg_qua_ka $
gov_eff_kai gov_eff_kao gov_eff_eq gov_eff_debt gov_eff_di reg_qua_kai reg_qua_kao reg_qua_eq reg_qua_debt $
reg_qua_di mm_flowg mmi_inflowg mmo_outflowg ca_gdp gov_debt_gdp age_dep_youth age_dep_old $
pop_grow ab_ca_gdp lca_ka lca_kai lca_kao lca_eq lca_debt lca_di reer_cpi tot ka_mean high_ka centile25 $
centile75
*
set spread = ka
set thresh = spread
linreg spread
# constant spread{1 2}
@UniqueValues(values=tvalues) thresh %regstart() %regend()
*
compute n=%rows(tvalues)
*
* These are the bottom and top of the permitted index values for the
* lower index, and the top of the permitted values for the upper index.
*
compute lstart=1
compute lend =n
*
set dtot = tot-tot{1}
set dreer_cpi = reer_cpi-reer_cpi{1}
set dca_gdp = ca_gdp-ca_gdp{1}
*
sweep
# dtot dreer_cpi dca_gdp
# constant dtot{1} dreer_cpi{1} dca_gdp{1}
compute loglbest=%logl
*
do lindex=lstart,lend
sweep(group=(thresh<=tvalues(lindex)))
# dtot dreer_cpi dca_gdp
# constant dtot{1} dreer_cpi{1} dca_gdp{1} spread{1}
if %logl>loglbest
compute lindexbest=lindex,loglbest=%logl
* dis lindexbest
end do lindex
disp "Best Break Value" tvalues(lindexbest)
*
compute lower=tvalues(lindexbest)
dec frml[int] switch
frml switch = 1+fix(ka>=lower)
*
* Estimate the model at the best breaks to get the covariance matrix.
*
sweep(group=switch(t))
# dtot dreer_cpi dca_gdp
# constant dtot{1} dreer_cpi{1} dca_gdp{1}
compute tvecmsigma=%sigma
*
diff tot / dtot
diff reer_cpi / dreer_cpi
diff ca_gdp / dca_gdp
*
system(model=basevecm)
variables dtot dreer_cpi dca_gdp
lags 1
det constant
end(system)
*
dec rect[frml] tvecfrml(3,2)
do i=1,2
estimate(smpl=(switch(t)==i))
frml(equation=%modeleqn(basevecm,1)) tvecfrml(1,i)
frml(equation=%modeleqn(basevecm,2)) tvecfrml(2,i)
frml(equation=%modeleqn(basevecm,3)) tvecfrml(3,i)
end do i
*
frml(identity) totid tot = tot{1}+dtot
frml(identity) reer_cpiid reer_cpi = reer_cpi{1}+dreer_cpi
frml(identity) ca_gdpid ca_gdp = ca_gdp{1}+dca_gdp
*
frml dtoteq dtot = tvecfrml(1,switch(t))
frml dreer_cpieq dreer_cpi = tvecfrml(2,switch(t))
frml dca_gdpeq dca_gdp = tvecfrml(3,switch(t))
*
group tvecm dtoteq dreer_cpieq dca_gdpeq totid reer_cpiid ca_gdpid
*
* Eventual forecast function, starting with 1995:1 data (largest value
* of spread).
*
*forecast(model=tvecm,from=2000:1,steps=10,results=eff)
*graph(footer=$
* "Eventual Forecast Function for Change in terms of trade, starting at 2000:1") 2
*# eff(1)
*# dtot 1//2000:1 1//2009:1
*graph(header=$
* "Eventual Forecast Function for Change in current account, starting at 2000:1") 2
*# eff(3)
*# dca_gdp 1//2000:1 1//2009:1
*graph(footer=$
* "Eventual Forecast Function for Change in real exchange rate, starting at 2000:1") 2
*# eff(2)
*# dreer_cpi 1//2000:1 1//2009:1
*
*graph(header=$
* "Eventual Forecast Function for current account, starting at 2000:1") 2
*# eff(6)
*# ca_gdp 2000:1 2009:1
*
* GIRF starting in 1//2000:1 for a one unit shock
* using the estimated covariance matrix. (2000:1 is in low ka regime)
*
compute ndraws=5000
compute baseentry=2000:1
compute nsteps =10
*
dec vect[series] fshocks(3)
dec vect[series] girf1(6) girf2(6) girf3(6)
dec vect[series] girf1_up(6) girf2_up(6) girf3_up(6)
dec vect[series] girf1_low(6) girf2_low(6) girf3_low(6)
dec series[vect] bishocks
dec vect ishocks
dec rect[series] sd_er(ndraws,6)
*
smpl baseentry baseentry+(nsteps-1)
*
do i=1,6
set girf3(i) = 0.0
set girf3_up(i) = 0.0
set girf3_low(i) = 0.0
end do i
*
do i=1,6
set girf2(i) = 0.0
set girf2_up(i) = 0.0
set girf2_low(i) = 0.0
end do i
*
do i=1,6
set girf1(i) = 0.0
set girf1_up(i) = 0.0
set girf1_low(i) = 0.0
end do i
*
*
compute fsigma=%psdfactor(tvecmsigma,||3,2,1||)
*
****
*girf for dca_gdp
*
do draw=1,ndraws
*random draw from 2000 to 2009 for each variable
gset bishocks = %ranmvnormal(fsigma)
set fshocks(1) = bishocks(t)(1)
set fshocks(2) = bishocks(t)(2)
set fshocks(3) = bishocks(t)(3)
*
*forecast(paths): make forecast using a series of shocks: "fshocks" in this case.
forecast(paths,model=tvecm,results=basesims)
# fshocks
*
*change initial period shock to the 3rd variable(dca_gdp) to one unit shock
compute ishock=%ran(fsigma(3,3))
*redo forcast after the change
compute fshocks(1)(baseentry)=bishocks(baseentry)(1)
compute fshocks(2)(baseentry)=bishocks(baseentry)(2)
compute fshocks(3)(baseentry)=bishocks(baseentry)(3)+ishock
forecast(paths,model=tvecm,results=sims)
# fshocks
do i=1,6
set girf3(i) = girf3(i)+(sims(i)-basesims(i))/ishock
set sd_er(draw,i) = (sims(i)-basesims(i))/ishock
end do i
end do draw
*
do i=1,6
set girf3(i) = girf3(i)/ndraws
end do i
*
* 84% upper limit, 16% lower limit
*
do i=1,6
do t=baseentry,baseentry+nsteps-1
do draw=1,ndraws
compute girf3_up(i)(t)=girf3_up(i)(t)+(sd_er(draw,i)(t)-girf3(i)(t))^2
end do draw
compute girf3_low(i)(t)=girf3(i)(t)-sqrt(girf3_up(i)(t)/(ndraws-1))
compute girf3_up(i)(t)=girf3(i)(t)+sqrt(girf3_up(i)(t)/(ndraws-1))
end do t
end do i
*
****
*girf for dreer_cpi
*
do draw=1,ndraws
*random draw from 2000 to 2009 for each variable
gset bishocks = %ranmvnormal(fsigma)
set fshocks(1) = bishocks(t)(1)
set fshocks(2) = bishocks(t)(2)
set fshocks(3) = bishocks(t)(3)
*
*forecast(paths): make forecast using a series of shocks: "fshocks" in this case.
forecast(paths,model=tvecm,results=basesims)
# fshocks
*
*change initial period shock to the 2nd variable to one sd shock
compute ishock=%ran(fsigma(2,2))
*
*redo forcast after the change
compute fshocks(1)(baseentry)=bishocks(baseentry)(1)
compute fshocks(2)(baseentry)=bishocks(baseentry)(2)+ishock
compute fshocks(3)(baseentry)=bishocks(baseentry)(3)
forecast(paths,model=tvecm,results=sims)
# fshocks
do i=1,6
set girf2(i) = girf2(i)+(sims(i)-basesims(i))/ishock
set sd_er(draw,i) = (sims(i)-basesims(i))/ishock
end do i
end do draw
*
do i=1,6
set girf2(i) = girf2(i)/ndraws
end do i
*
*
* 84% upper limit, 16% lower limit
*
do i=1,6
do t=baseentry,baseentry+nsteps-1
do draw=1,ndraws
compute girf2_up(i)(t)=girf2_up(i)(t)+(sd_er(draw,i)(t)-girf2(i)(t))^2
end do draw
compute girf2_low(i)(t)=girf2(i)(t)-sqrt(girf2_up(i)(t)/(ndraws-1))
compute girf2_up(i)(t)=girf2(i)(t)+sqrt(girf2_up(i)(t)/(ndraws-1))
end do t
end do i
*
****
*girf for dtot
*
do draw=1,ndraws
*random draw from 2000 to 2009 for each variable
gset bishocks = %ranmvnormal(fsigma)
set fshocks(1) = bishocks(t)(1)
set fshocks(2) = bishocks(t)(2)
set fshocks(3) = bishocks(t)(3)
*
*forecast(paths): make forecast using a series of shocks: "fshocks" in this case.
forecast(paths,model=tvecm,results=basesims)
# fshocks
*
*change initial period shock to the 1st variable(dca_gdp) to one unit shock
compute ishock=%ran(fsigma(1,1))
*
*redo forcast after the change
compute fshocks(1)(baseentry)=bishocks(baseentry)(1)+ishock
compute fshocks(2)(baseentry)=bishocks(baseentry)(2)
compute fshocks(3)(baseentry)=bishocks(baseentry)(3)
forecast(paths,model=tvecm,results=sims)
# fshocks
do i=1,6
set girf1(i) = girf1(i)+(sims(i)-basesims(i))/ishock
set sd_er(draw,i) = (sims(i)-basesims(i))/ishock
end do i
end do draw
*
do i=1,6
set girf1(i) = girf1(i)/ndraws
end do i
*
*
*
* 84% upper limit, 16% lower limit
*
do i=1,6
do t=baseentry,baseentry+nsteps-1
do draw=1,ndraws
compute girf1_up(i)(t)=girf1_up(i)(t)+(sd_er(draw,i)(t)-girf1(i)(t))^2
end do draw
compute girf1_low(i)(t)=girf1(i)(t)-sqrt(girf1_up(i)(t)/(ndraws-1))
compute girf1_up(i)(t)=girf1(i)(t)+sqrt(girf1_up(i)(t)/(ndraws-1))
end do t
end do i
*
spgraph(vfield=3,hfield=3,vlabel="response of",header="Low capital controls regime GIRF",$
footer="16% to 84% error bands")
graph(row=1,col=1,vlabel="change in terms of trade",header="one S.D. shock to terms of trade") 3
#girf1(1)
#girf1_up(1) * * 2
#girf1_low(1) * * 2
graph(row=2,col=1,vlabel="change in real exchange rate") 3
#girf1(2)
#girf1_up(2) * * 2
#girf1_low(2) * * 2
graph(row=3,col=1,vlabel="change in current account") 3
#girf1(3)
#girf1_up(3) * * 2
#girf1_low(3) * * 2
graph(row=1,col=2,header="one S.D. shock to real exchange rate") 3
#girf2(1)
#girf2_up(1) * * 2
#girf2_low(1) * * 2
graph(row=2,col=2) 3
#girf2(2)
#girf2_up(2) * * 2
#girf2_low(2) * * 2
graph(row=3,col=2) 3
#girf2(3)
#girf2_up(3) * * 2
#girf2_low(3) * * 2
graph(row=1,col=3,header="one S.D. shock to current account") 3
#girf3(1)
#girf3_up(1) * * 2
#girf3_low(1) * * 2
graph(row=2,col=3) 3
#girf3(2)
#girf3_up(2) * * 2
#girf3_low(2) * * 2
graph(row=3,col=3) 3
#girf3(3)
#girf3_up(3) * * 2
#girf3_low(3) * * 2
spgraph(done)