The RATS Software Forum

Code: Select all

*
open data el_ng_2014_june_10_2.xls
calendar(q) 1988:1
all 2013:1
data(format=xls, org=columns)
table
*******
*
@varlagselect(Lags=6,crit=hq)
# loil lcoal lnatgas lelhous
@mssysregression(states=2,switch=ch)
# loil lcoal lnatgas lelhous
# constant loil{1 to 1} lcoal{1 to 1} lnatgas{1 to 1} lelhous{1 to 1}
*
gset pt_t    = %zeros(nstates,1)
gset pt_t1   = %zeros(nstates,1)
gset psmooth = %zeros(nstates,1)
*
compute gstart=1989:1,gend=2013:1
@MSSysRegInitial gstart gend
compute p=%mspexpand(p)
gset MSRegime gstart gend = %ranbranch(||.5,.5||)
*
* Prior for transitions. Weak Dirichlet priors with preference for
* staying in a given regime.
*
dec vect[vect] gprior(nstates)
compute gprior(1)=||8.0,2.0||
compute gprior(2)=||2.0,8.0||
dec vect tcounts(nstates)
*
* Hierarchical prior for sigmas
* Uninformative prior on the common component.
*
* Priors for the relative Wisharts
*
dec vect nuprior(nstates)
compute nuprior=%fill(nstates,1,6.0)
*
dec vect[series] vresids(nvar)
dec vect[symm]   uu(nstates)
*
* Flat prior for the regressions.
*
dec symm hprior
dec vect bprior
compute hprior=%zeros(nvar*nreg,nvar*nreg)
compute bprior=%zeros(nvar*nreg,1)
*
compute nburn=1000, ndraws=10000
nonlin(parmset=allparms) betasys sigmav p
*
* For relabeling
*
dec vect voil(nstates)
dec vect[rect] tempbeta(nstates)
dec vect[symm] tempsigmav(nstates)
dec rect       tempp(nstates,nstates)
*
* Bookkeeping arrays
*
dec series[vect] bgibbs
gset bgibbs 1 ndraws = %parmspeek(allparms)
set regime1 gstart gend = 0.0
*
* Bookkeeping information for impulse responses
*
declare vect[rect] %%responses
dim %%responses(ndraws)
compute steps=15
*
infobox(action=define,lower=-nburn,upper=ndraws,progress) $
   "Gibbs Sampling"
do draw=-nburn,ndraws
   *
   * Draw sigma's given beta's and regimes
   * Compute the regime-specific residuals
   *
   @MSSysRegResids(regime=MSRegime) vresids gstart gend
   *
   * Compute the sums of squared residuals for each regime.
   *
   ewise uu(i)=%zeros(nvar,nvar)
   ewise tcounts(i)=0.0
   do time=gstart,gend
      compute uu(MSRegime(time))=uu(MSRegime(time))+$
          %outerxx(%xt(vresids,time))
      compute tcounts(MSRegime(time))=tcounts(MSRegime(time))+1
   end do time
   compute uucommon=%zeros(nvar,nvar)
   do k=1,nstates
      compute uucommon=uucommon+nuprior(k)*inv(sigmav(k))
   end do k
   *
   * Draw the common sigma given the regime-specific ones.
   *
   compute sigma=%ranwishartf(%decomp(inv(uucommon)),%sum(nuprior))
   *
   * Draw the regime-specific sigmas given the common one.
   *
   do k=1,nstates
      compute sigmav(k)=%ranwisharti(%decomp(inv(uu(k)+nuprior(k)*sigma)),$
          tcounts(k)+nuprior(k))
   end do k
   *
   * Draw beta's given sigma's and regimes
   *
   do i=1,nstates
      cmom(smpl=(MSRegime==i),model=MSSysRegModel) gstart gend
      compute betasys(i)=%reshape($
          %ranmvkroncmom(%cmom,inv(sigmav(i)),hprior,bprior),nreg,nvar)
   end do i
   *
   *
   *
   ewise voil(i)=sigmav(i)(3,3)
   compute swaps=%index(voil)
   if swaps(1)==2
      disp "Draw" draw "Executing swap"
   *
   * Relabel the beta's
   *
   ewise tempbeta(i)=betasys(i)
   ewise betasys(i)=tempbeta(swaps(i))
   *
   * Relabel the sigma's
   *
   ewise tempsigmav(i)=sigmav(i)
   ewise sigmav(i)=tempsigmav(swaps(i))
   *
   * Relabel the transitions
   *
   compute tempp=p
   ewise p(i,j)=tempp(swaps(i),swaps(j))
   *
   * Draw the regimes
   *
:skipbeta
   @MSSysRegFilter gstart gend
:redrawregimes
   @%mssample(start=gstart,end=gend) p pt_t pt_t1 MSRegime
   do i=1,nstates
      sstats gstart gend MSRegime==i>>tcounts(i)
   end do i
   if %minvalue(tcounts)<5 {
      disp "Draw" draw "Redrawing regimes with regime of size" $
         %minvalue(tcounts)
      goto redrawregimes
   }
   *
   * Draw p's
   *
   do j=1,nstates
      do i=1,nstates
         sstats(smpl=(MSRegime{1}==j.and.MSRegime==i)) gstart+1 gend $
            1>>tcounts(i)
      end do i
      compute pdraw=%randirichlet(tcounts+gprior(j))
      do i=1,nstates
         compute p(i,j)=pdraw(i)
      end do i
   end do j
   infobox(current=draw)
   if draw>0 {
      *
      * Do the bookkeeping
      *
      set regime1 gstart gend = regime1+(MSRegime==1)
      compute bgibbs(draw)=%parmspeek(allparms)
      *
      * Save IRF's for regime 1. Choleski factor standardized to unit
      * shock is used.
      *
      @MSSysRegSetModel(regime=1)
      compute factor=%decomp(sigmav(1)),factor=factor*inv(%diag(%xdiag(factor)))
      impulse(noprint,model=MSSysRegModel,results=impulses1,steps=steps,factor=factor)
      *
      * Save IRF's for regime 2
      *
      @MSSysRegSetModel(regime=2)
      compute factor=%decomp(sigmav(2)),factor=factor*inv(%diag(%xdiag(factor)))
      impulse(noprint,model=MSSysRegModel,results=impulses2,steps=steps,factor=factor)
      *
      * This will interleave the shocks from the two regimes.
      *
      dim %%responses(draw)(%rows(impulses1)*%cols(impulses1)*2,steps)
      ewise %%responses(draw)(i,j)=ix=%vec(%xt(impulses1,j)~~%xt(impulses2,j)),ix(i)
   }
end do draw
infobox(action=remove)
@mcmcpostproc(ndraws=ndraws,mean=bmeans,stderrs=bstderrs) bgibbs
*
report(action=define)
report(atrow=1,atcol=1,fillby=cols) %parmslabels(allparms)
report(atrow=1,atcol=2,fillby=cols) bmeans
report(atrow=1,atcol=3,fillby=cols) bstderrs
report(action=format,picture="*.###")
report(action=show)
*
set regime1 gstart gend = regime1/ndraws
graph(header="MCMC Probability of  Regime 1")
# regime1
*
@MCProcessIRF(model=MSSysRegModel,lower=lower,upper=upper,irf=irf,percentile=||.025,.975||)

table
set irf4_1 = IRF(4,1)
set lower4_1 = LOWER(4,1)
set upper4_1 = UPPER(4,1)
*

set irf4_2 = IRF(4,2)
set lower4_2 = LOWER(4,2)
set upper4_2 = UPPER(4,2)
*
set irf4_3 = IRF(4,3)
set lower4_3 = LOWER(4,3)
set upper4_3 = UPPER(4,3)

*
set irf4_4 = IRF(4,4)
set lower4_4 = LOWER(4,4)
set upper4_4 = UPPER(4,4)
*
set irf4_4 = IRF(4,4)
set lower4_4 = LOWER(4,4)
set upper4_4 = UPPER(4,4)

*
set irf4_6 = IRF(4,6)
set lower4_6 = LOWER(4,6)
set upper4_6 = UPPER(4,6)

*
set irf4_7 = IRF(4,7)
set lower4_7 = LOWER(4,7)
set upper4_7 = UPPER(4,7)

*
set irf4_8 = IRF(4,8)
set lower4_8 = LOWER(4,8)
set upper4_8 = UPPER(4,8)

*******************************************
set regime2 = 1- regime1
spgraph(hfields=1,vfields=2)
graph(header="Regime 1", style=BAR, max=1.0)
# regime1
graph(header="Regime 2", style=BAR,max=1.0)
# regime2
spgraph(done)
print / regime1 regime2
*****************Linear VAR Model**************************
system(model=modelin)
variables loil lcoal lnatgas lelhous
lags 1 to 2
det constant
end(system)
estimate(resids=varresids, noprint)
errors(model=modelin,window="Variance Decomposition(Linear VAR)", steps=15)
impulse(model=modelin,window="Impulse Responses",steps=40)
source mcvardodrawsunit.src
@MCVARDoDrawsunit(model=modelin, steps=40)
@MCProcessIRF(model=modelin,irf=irflin,center=mean, lower=lowerlin, upper = upperlin,percentile=||.025,.975||)

* saving responses and error bands
* responses to ip under different regimes
* linear model
set ip1_lin 1 15 = irflin(1,1)
set ip2_lin 1 15 = irflin(2,1)
set ip3_lin 1 15 = irflin(3,1)
set ip4_lin 1 15 = irflin(4,1)
print / ip5_lin
* saving upper and lower error bands
set ip1_lin_up 1 15 UPPERLIN(1,1)
set ip2_lin_up 1 15 UPPERLIN(2,1)
set ip3_lin_up 1 15 = UPPERLIN(3,1)
set ip4_lin_up 1 15 UPPERLIN(4,1)


set ip1_lin_low 1 15 = LOWERLIN(1,1)
set ip2_lin_low 1 15 LOWERLIN(2,1)
set ip3_lin_low 1 15 = LOWERLIN(3,1)
set ip4_lin_low 1 15 LOWERLIN(4,1)


* intrate linear model
set int1_lin 1 15 irflin(1,2)
set int2_lin 1 15 = irflin(2,2)
set int3_lin 1 15 = irflin(3,2)
set int4_lin 1 15 = irflin(4,2)

* saving upper and lower error bands for cpi
set int1_lin_up 1 15 = UPPERLIN(1,2)
set int2_lin_up 1 15 = UPPERLIN(2,2)
set int3_lin_up 1 15 = UPPERLIN(3,2)
set int4_lin_up 1 15 = UPPERLIN(4,2)


set int1_lin_low 1 15 = LOWERLIN(1,2)
set int2_lin_low 1 15 = LOWERLIN(2,2)
set int3_lin_low 1 15 = LOWERLIN(3,2)
set int4_lin_low 1 15 = LOWERLIN(4,2)

*
* linear model M1
set m1_lin 1 15 = irflin(1,3)
set m2_lin 1 15 = irflin(2,3)
set m3_lin 1 15 = irflin(3,3)
set m4_lin 1 15 = irflin(4,3)



* saving upper and lower error bands for int
set m1_lin_up 1 15 = UPPERLIN(1,3)
set m2_lin_up 1 15 = UPPERLIN(2,3)
set m3_lin_up 1 15 = UPPERLIN(3,3)
set m4_lin_up 1 15 = UPPERLIN(4,3)



*
set m1_lin_low 1 15 = LOWERLIN(1,3)
set m2_lin_low 1 15 = LOWERLIN(2,3)
set m3_lin_low 1 15 = LOWERLIN(3,3)
set m4_lin_low 1 15 = LOWERLIN(4,3)


*
* responses to cp under different regimes
* linear model
set cp1_lin 1 15 = irflin(1,4)
set cp2_lin 1 15 = irflin(2,4)
set cp3_lin 1 15 = irflin(3,4)
set cp4_lin 1 15 = irflin(4,4)

* saving upper and lower error bands for cpi
set cp1_lin_up 1 15 = UPPERLIN(1,4)
set cp2_lin_up 1 15 = UPPERLIN(2,4)
set cp3_lin_up 1 15 = UPPERLIN(3,4)
set cp4_lin_up 1 15 = UPPERLIN(4,4)

*
set cp1_lin_low 1 15 = LOWERLIN(1,4)
set cp2_lin_low 1 15 = LOWERLIN(2,4)
set cp3_lin_low 1 15 = LOWERLIN(3,4)
set cp4_lin_low 1 15 = LOWERLIN(4,4)


*
GRPARM(FONT="Arial") SUBHEADER 32
GRPARM(FONT="Arial") HEADER 34
GRPARM(FONT="Arial Narrow") AXISLABELS 24
GRPARM(FONT="Arial") AXISLABELS 28


spgraph(hfields=3,vfields=4)
graph(nodate, header ="Responses to Oil Price Shocks: Linear VAR", row=1, col=1, min=-.4, max=1.2)  3
# ip4_lin  * 40
# ip4_lin_up * 40 2
# ip4_lin_low * 40 2
graph(nodate, header ="Responses to Oil Price Shocks: Regime 1", row=1, col=2, min=-.4, max=1.2) 3
# irf4_1 * 40
# lower4_1 * 40  2
# upper4_1 * 40  2
*
graph(nodate,header ="Responses to Oil Price Shocks: Regime 2",row=1, col=3, min=-.4, max=1.2) 3
# irf4_2 * 40
# lower4_2 * 40 2
# upper4_2 * 40 2

*interest rates
graph(nodate, header ="Responses to Coal Price Shocks: Linear VAR",row=2, col=1, min=-.075, max=.075) 3
# int4_lin  * 40
# int4_lin_up * 40 2
# int4_lin_low * 40 2

*
graph(nodate,header ="Responses to Coal Price Shocks: Regime 1",row=2, col=2) 3
# irf4_3 * 40
# lower4_3 * 40 2
# upper4_3 * 40 2
*
graph(nodate,header ="Responses to Coal Price Shocks: Regime 2" , row=2, col=3, min=-.075, max=.075) 3
# irf4_4 * 40
# lower4_4 * 40 2
# upper4_4 * 40 2

* m1 shocks
graph(nodate, header ="Responses to Natural Gas Price Shocks: Linear VAR",row=3, col=1, min=-.5, max=1.0) 3
# m4_lin  * 40
# m4_lin_up  * 40 2
# m4_lin_low * 40 2
*
graph(nodate,header ="Responses to Natural Gas Price Shocks: Regime 1",row=3, col=2, min=-.5, max=1.0) 3
# irf4_5  * 40
# lower4_5 * 40 2
# upper4_5 * 40 2
**
graph(nodate,header ="Responses to Natural Gas Price Shocks: Regime 2", row=3, col=3, min=-.5, max=1.0) 3
# irf4_6 * 40
# lower4_6 * 40 2
# upper4_6 * 40 2
*cpi shocks
graph(nodate, header="Responses to Electricity Price Shocks: Linear VAR", row=4, col=1, min=-1.5, max=1.5) 3
# cp4_lin * 40
# cp4_lin_up * 40 2
# cp4_lin_low  * 40 2
*
graph(nodate,header ="Responses to Electricity Price Shocks: Regime 1", row=4, col=2, min=-1.5, max=1.5) 3
# irf4_7 * 40
# lower4_7 * 40 2
# upper4_7 * 40 2
*
graph(nodate,header ="Responses to Electricity Price Shocks: Regime 2", row=4, col=3, min=-1.5, max=1.5) 3
# irf4_8 * 40
# lower4_8 * 40 2
# upper4_8 * 40 2
spgraph(done)

Code: Select all

ewise voil(i)=sigmav(i)(3,3)
compute swaps=%index(voil)
if swaps(1)==2
disp "Draw" draw "Executing swap"