Ehrmann-Ellison-Valla(2003) Regime-dependent IRF's

[ege_man]

Dear Tom

I have no outlier in the data since I am using level of the variables in the estimation. Here is the code modified for my study. The data is also attached. I would really appreciate if you tell me what is wrong with the estimates.

Regards

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)

[TomDoan]

Why do you think you have no outliers? Have you run just a straight VAR? Your oil series has a couple of almost 5 standard error residuals, and the others have 3's and 4's. If you allow variances to switch, it's going to try to fix those.

[ege_man]

Yes I run a linear VAR model if you look at the end of the code I obtained positive and significant responses but standard error bands are large. So do you mean that MS-VAR estimation is not applicable for this data? Do you have any suggestion to model regime change, because I would like to capture the impact of privatizations in the electricity market.

[TomDoan]

Why would you think a Markov-switching model would be appropriate to study that? In the MS model, the regime is unobservable, and can switch at any time from one to the other. Wouldn't privatization take place at known times?

[allister]

Hi Tom this may be a very dumb question but i will ask nonetheless. Is it possible to estimate this model with two (2) variables and if so what are the parts of the code that i will have to amend?


Thanks in advance for your guidance.

[TomDoan]

Change this (which you would have to do to match your variables anyway)

@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}

and any graphics.

[aloprofundo]

Dear Tom,

I have estimated the model with EM method, now I need to compute the confidence bands for the IRF. It seems that it is possible to obtain the bands using MCMC method, as the following paper has done:
http://web.up.ac.za/sitefiles/file/40/6 ... 012_22.pdf

Would you please provide any clue how to do it?
Thanks!!

I'm not sure where you got the code for using flipper, but you're multiplying it by a matrix that isn't used in the impulse responses. What you want is something like:

compute flipper=||1.0,1.0,-1.0||
compute factor=%decomp(sigmav(2)),factor=factor*inv(%diag(%xdiag(factor)))
compute factor=factor*%diag(flipper)

As written, this will make the impact of shock 3 on the interest rate negative, and the impact of shock 1 on hpp and shock 2 on gdpp positive.

If you want to do confidence bands, use the MCMC version.

[TomDoan]

Dear Tom,

I have estimated the model with EM method, now I need to compute the confidence bands for the IRF. It seems that it is possible to obtain the bands using MCMC method, as the following paper has done:
http://web.up.ac.za/sitefiles/file/40/6 ... 012_22.pdf

Would you please provide any clue how to do it?
Thanks!!

I'm not sure where you got the code for using flipper, but you're multiplying it by a matrix that isn't used in the impulse responses. What you want is something like:

compute flipper=||1.0,1.0,-1.0||
compute factor=%decomp(sigmav(2)),factor=factor*inv(%diag(%xdiag(factor)))
compute factor=factor*%diag(flipper)

As written, this will make the impact of shock 3 on the interest rate negative, and the impact of shock 1 on hpp and shock 2 on gdpp positive.

If you want to do confidence bands, use the MCMC version.

They use some form of "bootstrapping" (somewhat ill-defined). As with the advice that you quoted, if you want confidence bands, use the MCMC version.

[aloprofundo]

Hi Tom

I realized that the EM version is more suitable for my research, because it allows me to do the specification test while the Bayesian method doesn't (correct me if I wrong). Moreover, I am not really a Bayesian person so it would take me quite a while to understand its background if I choose the MCMC version. I am simply wondering if you could provide some hints on using the bootstrapping technique to construct the impulse response's confidence interval, as the following discussion paper (Ehrmann et al., 2001) suggests:

http://www.suomenpankki.fi/en/julkaisut ... s/0111.pdf (see p.14 )

I really appreciate your help!

[TomDoan]

That's more complicated than doing the MCMC.

[aloprofundo]

Hi Tom, sorry to bother again, but I wonder if it would be possible to use EM to estimate MS-VAR model and apply MCMC to construct impulse responses' confidence intervals, if the answer is yes, could you provide a hint (which part of the code to look at) how to do so?

[TomDoan]

Hi Tom, sorry to bother again, but I wonder if it would be possible to use EM to estimate MS-VAR model and apply MCMC to construct impulse responses' confidence intervals, if the answer is yes, could you provide a hint (which part of the code to look at) how to do so?

No. I'm not sure I understand the reluctance to do MCMC. As modern Bayesian methods go, this is fairly straightforward, since it's all sampling from known distributions. MCMC is really the only good way to attack Markov switching models with differing variances because the likelihood is unbounded.

[zy761]

Dear Tom,

If I understand right, in this program you show the oil shock (3rd variable) to all the 3 variables in the system. I notice the following in your program:

Code: Select all

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

Do I need to change sigmav(i)(3,3) with other numbers if I want to estimate shock from a different variable (not the 3rd one)?

As I am not quite clear the role of this part, any clarification and help would be very appreciated. Thank you!

[TomDoan]

That's being used to "label" the regimes based upon the variance of the oil shock (in this case). The biggest problem with MS models with so many parameters free is that the regimes can switch labels.

[PTillmann-436]

Dear Tom,

I need to estimate regime-dependent IRFs for shocks to an exogenous variable, i.e. an MS-VAR-X, such as here:

https://estima.com/forum/viewtopic.php?f=4&t=1154#p4085

Is there a ways to incluce a placeholder equation for the exogenous variable such as in the example? Do you have another idea of how to handle it?

Thank you so much

Peter