Identifying VARs with sign restrictions

[TomDoan]

There's an article in the forthcoming newsletter on this very topic. See

http://www.estima.com/newslett/March2009RATSLetter.pdf

[fabio fornari]

thanks a lot Tom.
will see how it works with three shocks before moving on ... needs really lots of draws to get acceptances ...
will let you know
thanks f

[mmk]

Hi

Did anyone try to extend Uhilg_3 file for the paper published (2005) with penalty function to more than 6 variables? can anyone tell me how the penalty is going to change in this case?

Thanks for help

[TomDoan]

Hi

Did anyone try to extend Uhilg_3 file for the paper published (2005) with penalty function to more than 6 variables? can anyone tell me how the penalty is going to change in this case?

Thanks for help

This uses an alternative and more flexible parameterization of the unit sphere than was used by Uhlig. (Note that it's really the parameterization that needs to change; the calculation of the penalty function is similar regardless). I didn't use it in the original uhlig3.prg because the %stereo function was added with version 6.35 of RATS. There are quite a few ways to parameterize the unit sphere. The stereoscopic projection used in the %stereo function takes the N-1 vector x and maps it to the N-vector formed by

y(i)=2x(i)/(||x||^2+1)
y(N)=(||x||^2-1)/(||x||^2+1)

All such mappings have at least one minor flaw. In this case, the only way to make y(N)=0 is to make the x's very large. (Uhlig's mapping has the problem of non-uniqueness, as well as being rather difficult to generalize easily).

Code: Select all

*
* Replication File for Uhlig(2005), "What are the effects of monetary policy on
* output? Results from an agnostic identification procedure," Journal of Monetary
* Economics, 52, pp 381-419. Penalty function approach.
*
open data uhligdata.xls
calendar(m) 1965
data(format=xls,org=columns) 1965:1 2003:12 gdpc1 gdpdef cprindex totresns bognonbr fedfunds
*
set gdpc1    = log(gdpc1)*100.0
set gdpdef   = log(gdpdef)*100.0
set cprindex = log(cprindex)*100.0
set totresns = log(totresns)*100.0
set bognonbr = log(bognonbr)*100.0
*
system(model=varmodel)
variables gdpc1 gdpdef cprindex fedfunds bognonbr totresns
lags 1 to 12
end(system)
estimate(noprint)
*
dec vect[strings] vl(6)
compute vl=||"real GDP","GDP price defl","Comm. Price Ind.",$
             "Fed Funds Rate","Nonborr. Reserv.","Total Reserves"||
*
* Get the scale factors to be used for impulse responses
*
dec vect scales(6)
dec real func
compute fill=0
dofor i = gdpc1 gdpdef cprindex fedfunds bognonbr totresns
   diff i / diff1
   stats(noprint) diff1
   compute fill=fill+1,scales(fill)=sqrt(%variance)
end dofor i
*
* ndraws is the desired number of draws from the posterior of the VAR
* nvar is the number of variables
* nstep is the number of IRF steps to compute
* KMAX is the "K" value for the number of steps constrained
*
*
compute ndraws=1000
compute nvar=6
compute nstep=60
compute KMAX=5
declare vect[rect] goodresp(ndraws)
declare vector ik a(nvar)
*
dec vect g(nvar-1)
compute g=%fill(nvar-1,1,1.0)
nonlin g
*
* This is the standard setup for MC integration of an OLS VAR
*
compute sxx    =%decomp(%xx)
compute svt    =%decomp(inv(%nobs*%sigma))
compute betaols=%modelgetcoeffs(varmodel)
compute ncoef  =%rows(sxx)
compute wishdof=%nobs-ncoef
dec rect ranc(ncoef,nvar)
*
infobox(action=define,progress,lower=1,upper=ndraws) "Monte Carlo Integration"
*
* Plan for double the number of draws to allow for the ones which we
* reject.
*
compute accept=0
do draws=1,ndraws*2
   *
   * Make a draw from the posterior for the VAR and compute its impulse
   * responses.
   *
   compute sigmad  =%ranwisharti(svt,wishdof)
   compute swish   =%decomp(sigmad)
   compute ranc    =%ran(1.0)
   compute betau   =sxx*ranc*tr(swish)
   compute betadraw=betaols+betau
   compute %modelsetcoeffs(varmodel,betadraw)
   impulse(noprint,model=varmodel,decomp=swish,results=impulses,steps=nstep)
   *
   * Minimize the penalty function, starting from the last set of minimizers
   *
   find(noprint) min func
      compute a=%stereo(g)
      compute func=0.0
      do k=1,KMAX+1
         compute ik=(%xt(impulses,k)*a)./scales
         compute func=func+%if(ik(4)>0,-ik(4),-100*ik(4))+$
                           %if(ik(3)<0,ik(3),100*ik(3))+$
                           %if(ik(2)<0,ik(2),100*ik(2))+$
                           %if(ik(5)<0,ik(5),100*ik(5))
      end do k
   end find
   compute testfunc=func,testbeta=%beta
   *
   * Try the minimization again, starting from a standard set of values
   *
   compute g=%fill(nvar-1,1,1.0)
   find(noprint) min func
      compute a=%stereo(g)
      compute func=0.0
      do k=1,KMAX+1
         compute ik=(%xt(impulses,k)*a)./scales
         compute func=func+%if(ik(4)>0,-ik(4),-100*ik(4))+$
                           %if(ik(3)<0,ik(3),100*ik(3))+$
                           %if(ik(2)<0,ik(2),100*ik(2))+$
                           %if(ik(5)<0,ik(5),100*ik(5))
      end do k
   end find
   *
   * If the two estimates don't match, reject the draw. If they do, copy out the
   * impulse responses.
   *
   if abs(testfunc-func)<.01 {
      compute accept=accept+1
      dim goodresp(accept)(nstep,nvar)
      ewise goodresp(accept)(i,j)=ik=%xt(impulses,i)*a,ik(j)
   }
   *
   * If we've hit out desired number of accepted draws, break
   *
   if accept>=ndraws
      break
   infobox(current=accept)
end do draws
infobox(action=remove)
*
* Post-processing. Graph the mean of the responses along with the 16% and 84%-iles
*
clear upper lower resp
*
spgraph(vfields=3,hfields=2,hlabel="Figure 14. Impulse Responses with Penalty Function Approach")
do i=1,nvar
   compute minlower=maxupper=0.0
   smpl 1 accept
   do k=1,nstep
      set work = goodresp(t)(k,i)
      compute frac=%fractiles(work,||.16,.84||)
      compute lower(k)=frac(1)
      compute upper(k)=frac(2)
      compute resp(k)=%avg(work)
   end do k
*
   smpl 1 nstep
   graph(ticks,number=0,picture="##.##",header="Impulse Responses for "+vl(i)) 3
   # resp
   # upper / 2
   # lower / 2
end do i
*
spgraph(done)


[mmk]

Many thanks for providing this flexible way

[mmk]

Hi

I tried to run the new uhilg_3 using the parameterization you suggested for the penalty function. It produces graphs that differ from the analysis of uhilg in his original paper. GDP deflator as well as PPI and nonborrowered reserves should be negative by construction during the first 6 months even they still take those ngative valuse soforth. Can you explain to me why they are different?

thanks for help

[lhlee0506]

Hi~Tom,

Would please you kindly instruct me how to collect the value of the initial demand/monetary shock and how to show the variance of decomposition?

[TomDoan]

Hi

I tried to run the new uhilg_3 using the parameterization you suggested for the penalty function. It produces graphs that differ from the analysis of uhilg in his original paper. GDP deflator as well as PPI and nonborrowered reserves should be negative by construction during the first 6 months even they still take those ngative valuse soforth. Can you explain to me why they are different?

thanks for help

Thanks for pointing that out. The code posted above has been fixed. In patching the code over, I forgot the new nonlin instruction below.

Code: Select all

dec vect g(nvar-1)
compute g=%fill(nvar-1,1,1.0)
nonlin g

[lhlee0506]

Hi~

I tried to put the variance decomposition in the previous code with 2 shocks you posted few weeks ago.
and I got message as following,

## SX22. Expected Type REAL, Got MATRIX(REAL) Instead
>>>>%xt(impulses,t).^2<<<<

Would you kindly instruct me how to fix it?

Thank you very much!

[TomDoan]

You need a

Code: Select all

DEC SERIES[RECT] irfsquared

instruction. Otherwise GSET doesn't know what type it's expecting and defaults to just a series of reals.

[lhlee0506]

Hello,

I tried to follow up your instruction in the latest RATS letter which demonstrateds the impact responses on variables 1 and 2 to be equal. But I failed. Perhaps I misunderstood your instruction.
If you could give me some guidances, I will be very appreciated.

Code: Select all

*
*  Replication File for Uhlig (2005), "What are the effects of monetary policy on output?
*  Results from an agnostic identification procedure." Journal of Monetary Economics, 52, pp
*  381-419. Pure sign restriction approach
*
open data levelsa1.xls
calendar 1970 1 4
data(format=xls,org=columns) 1970:01 2008:03 lrgspen1 lrntax lrgdp pgdpch rint lrpcon lrpinv
*

system(model=varmodel)
variables lrgspen1 lrntax lrgdp pgdpch rint lrpcon
lags 1 to 4
end(system)
estimate(noprint)
*
dec vect[strings] vl(6)
compute vl=||'Govt spending','net tax','real gdp','GDP deflator','real interest rate','Private onsumption'||
*
* n1 is the number of draws from the posterior of the VAR
* n2 is the number of draws from the unit sphere for each draw for the VAR
* nvar is the number of variables
* nstep is the number of IRF steps to compute
* KMAX is the "K" value for the number of steps constrained
*
compute n1=200
compute n2=200
compute nvar=6
compute nstep=60
compute KMAX=5
*
* This is the standard setup for MC integration of an OLS VAR
*
compute sxx    =%decomp(%xx)
compute svt    =%decomp(inv(%nobs*%sigma))
compute betaols=%modelgetcoeffs(varmodel)
compute ncoef  =%rows(sxx)
compute wishdof=%nobs-ncoef
dec rect ranc(ncoef,nvar)
*
* Most draws are going to get rejected. We allow for up to 1000
* good ones. The variable accept will count the number of accepted
* draws. GOODRESP will be a RECT(nsteps,nvar) at each accepted
* draw.
*
declare vect[rect] goodresp(1000)
declare vector ik a(nvar)
*
compute accept=0
infobox(action=define,progress,lower=1,upper=n1) 'Monte Carlo Integration'
do draws=1,n1
   *
   * Make a draw from the posterior for the VAR and compute its impulse
   * responses.
   *
   compute sigmad  =%ranwisharti(svt,wishdof)
   compute swish   =%decomp(sigmad)
   compute ranc    =%ran(1.0)
   compute betau   =sxx*ranc*tr(swish)
   compute betadraw=betaols+betau
   compute %modelsetcoeffs(varmodel,betadraw)
   impulse(noprint,model=varmodel,decomp=swish,results=impulses,steps=nstep)
   *
   * Do the subdraws over the unit sphere. These give the weights on the
   * orthogonal components.
   *
   do subdraws=1,n2
      compute rperp=$%perp(||1.0,-1.0,0.0,0.0,0.0,0.0||*swish)
      compute a=rperp*%ranmat(5,1),a=a/sqrt(%normsqr(a))
      *
      * Check that the responses have the correct signs for steps 1 to KMAX+1
      * (+1 because in the paper, the steps used 0-based subscripts, rather than
      * the 1 based subscripts used by RATS).
      *
      do k=1,KMAX+1
         compute ik=%xt(impulses,k)*a
         if ik(1)<0
         branch 105
      end do k
      *
      * This is an accepted draw. Copy the information out. If we have our 1000
      * good ones, drop out of the overall loop.
      *
      compute accept=accept+1
      dim goodresp(accept)(nstep,nvar)
      ewise goodresp(accept)(i,j)=ik=%xt(impulses,i)*a,ik(j)
      if accept>=1000
         break
   :105
   end do subdraws
   if accept>=1000
      break
   infobox(current=draws)
end do draws
infobox(action=remove)
*
* Post-processing. Graph the mean of the responses along with the 16% and 84%-iles
*
clear upper lower resp
*
spgraph(vfields=3,hfields=2,hlabel='Figure 6. Impulse Responses with Pure-Sign Approach',subhea='monetary Shock')
do i=1,nvar
   compute minlower=maxupper=0.0
   smpl 1 accept
   do k=1,nstep
      set work = goodresp(t)(k,i)
      compute frac=%fractiles(work,||.16,.84||)
      compute lower(k)=frac(1)
      compute upper(k)=frac(2)
      compute resp(k)=%avg(work)
   end do k
*
   smpl 1 nstep
   graph(ticks,number=0,picture='##.##',header='Impulse Responses for '+vl(i)) 3
   # resp
   # upper / 2
   # lower / 2
end do i
*
spgraph(done)















[TomDoan]

Hello,

I tried to follow up your instruction in the latest RATS letter which demonstrateds the impact responses on variables 1 and 2 to be equal. But I failed. Perhaps I misunderstood your instruction.
If you could give me some guidances, I will be very appreciated.

Code: Select all

      compute rperp=$%perp(||1.0,-1.0,0.0,0.0,0.0,0.0||*swish)
 

The line I isolated has been garbled. (Probably didn't fit in one column in the newsletter). Get rid of the $ and it should be fine.

[lhlee0506]

Hi~

Thank you for your instruction.
After removing $, the code seems to wotk.
But I am still confused about how to set up the sign restriction after putting variable1,2 to be equal.
This modification seems only make variable 1,2 to be equal for one period, how to make it last for four periods?
Hope to have your kind responses.
Thanks a lot!

Code: Select all

*
*  Replication File for Uhlig (2005), "What are the effects of monetary policy on output?
*  Results from an agnostic identification procedure." Journal of Monetary Economics, 52, pp
*  381-419. Pure sign restriction approach
*
open data levelsa1.xls
calendar 1970 1 4
data(format=xls,org=columns) 1970:01 2008:03 lrgspen1 lrntax lrgdp pgdpch rint lrpcon lrpinv
*

system(model=varmodel)
variables lrgspen1 lrntax lrgdp pgdpch rint lrpcon
lags 1 to 4
end(system)
estimate(noprint)
*
dec vect[strings] vl(6)
compute vl=||'Govt spending','net tax','real gdp','GDP deflator','real interest rate','Private onsumption'||
*
* n1 is the number of draws from the posterior of the VAR
* n2 is the number of draws from the unit sphere for each draw for the VAR
* nvar is the number of variables
* nstep is the number of IRF steps to compute
* KMAX is the "K" value for the number of steps constrained
*
compute n1=200
compute n2=200
compute nvar=6
compute nstep=60
compute KMAX=5
*
* This is the standard setup for MC integration of an OLS VAR
*
compute sxx    =%decomp(%xx)
compute svt    =%decomp(inv(%nobs*%sigma))
compute betaols=%modelgetcoeffs(varmodel)
compute ncoef  =%rows(sxx)
compute wishdof=%nobs-ncoef
dec rect ranc(ncoef,nvar)
*
* Most draws are going to get rejected. We allow for up to 1000
* good ones. The variable accept will count the number of accepted
* draws. GOODRESP will be a RECT(nsteps,nvar) at each accepted
* draw.
*
declare vect[rect] goodresp(1000)
declare vector ik a(nvar)
*
compute accept=0
infobox(action=define,progress,lower=1,upper=n1) 'Monte Carlo Integration'
do draws=1,n1
   *
   * Make a draw from the posterior for the VAR and compute its impulse
   * responses.
   *
   compute sigmad  =%ranwisharti(svt,wishdof)
   compute swish   =%decomp(sigmad)
   compute ranc    =%ran(1.0)
   compute betau   =sxx*ranc*tr(swish)
   compute betadraw=betaols+betau
   compute %modelsetcoeffs(varmodel,betadraw)
   impulse(noprint,model=varmodel,decomp=swish,results=impulses,steps=nstep)
   *
   * Do the subdraws over the unit sphere. These give the weights on the
   * orthogonal components.
   *
   do subdraws=1,n2
      compute rperp=%perp(||1.0,-1.0,0.0,0.0,0.0,0.0||*swish)
      compute a=rperp*%ranmat(5,1),a=a/sqrt(%normsqr(a))
      *
      * Check that the responses have the correct signs for steps 1 to KMAX+1
      * (+1 because in the paper, the steps used 0-based subscripts, rather than
      * the 1 based subscripts used by RATS).
      *
      do k=1,KMAX+1
         compute ik=%xt(impulses,k)*a
         if ik(1)<0
         branch 105
      end do k
      *
      * This is an accepted draw. Copy the information out. If we have our 1000
      * good ones, drop out of the overall loop.
      *
      compute accept=accept+1
      dim goodresp(accept)(nstep,nvar)
      ewise goodresp(accept)(i,j)=ik=%xt(impulses,i)*a,ik(j)
      if accept>=1000
         break
   :105
   end do subdraws
   if accept>=1000
      break
   infobox(current=draws)
end do draws
infobox(action=remove)
*
* Post-processing. Graph the mean of the responses along with the 16% and 84%-iles
*
clear upper lower resp
*
spgraph(vfields=3,hfields=2,hlabel='Figure 6. Impulse Responses with Pure-Sign Approach',subhea='monetary Shock')
do i=1,nvar
   compute minlower=maxupper=0.0
   smpl 1 accept
   do k=1,nstep
      set work = goodresp(t)(k,i)
      compute frac=%fractiles(work,||.16,.84||)
      compute lower(k)=frac(1)
      compute upper(k)=frac(2)
      compute resp(k)=%avg(work)
   end do k
*
   smpl 1 nstep
   graph(ticks,number=0,picture='##.##',header='Impulse Responses for '+vl(i)) 3
   # resp
   # upper / 2
   # lower / 2
end do i
*
spgraph(done)















[TL]

Dear Tom,

I use code for two shocks from the forum. In addition to the impulse responses, I would like to see how much of the variation is accounted for by the first shock and by the second shock (variance decomposition).

I am not quite sure whether fraction explained (variance decomposition) is correct. Could you please give me some suggestions and feedbacks whether my code is correct?

Thank you very much for your help and kindness.

Tim

Code: Select all

    OPEN DATA uhligdata.xls
    CALENDAR 1965 1 12
    compute missc=1.0e+32
    data(format=xlS,org=columns) 1965:1 2003:12 gdpc1 gdpdef cprindex totresns bognonbr fedfunds

    set gdpc1    = log(gdpc1)*100.0
    set gdpdef   = log(gdpdef)*100.0
    set cprindex = log(cprindex)*100.0
    set totresns = log(totresns)*100.0
    set bognonbr = log(bognonbr)*100.0
    *
    system(model=varmodel)
    variables gdpc1 gdpdef cprindex fedfunds bognonbr totresns
    lags 1 to 12
    end(system)
    estimate(noprint,resid=resids)

    dec vect[strings] vl(6)
    compute vl=||'real GDP','GDP price defl','Comm. Price Ind.',$
         'Fed Funds Rate','Nonborr. Reserv.','Total Reserves'||

    compute n1=4000
    compute n2=4000
    compute nvar=6
    compute nstep=60
    compute KMAX=5

    * This is the standard setup for MC integration of an OLS VAR
    *
    dec symm s(6,6)
    dec vect v1(6)      ;* For the unit vector on the 1st draw
    dec vect v2(5)      ;* For the unit vector on the 2nd draw
    dec vect v(6)     ;* Working impulse vector

    compute sxx    =%decomp(%xx)
    compute svt    =%decomp(inv(%nobs*%sigma))
    compute betaols=%modelgetcoeffs(varmodel)
    compute ncoef  =%rows(sxx)
    compute wishdof=%nobs-ncoef
    dec rect ranc(ncoef,nvar)
    *
    * Most draws are going to get rejected. We allow for up to 1000
    * good ones. The variable accept will count the number of accepted
    * draws. GOODRESP will be a RECT(nsteps,nvar) at each accepted
    * draw.
    *
    declare vect[rect] goodresp(1000) goodfevd(1000)
    declare vect[rect] goodrespa(1000) goodfevda(1000)

    declare vector ik a(nvar) ones(nvar)
    declare series[rect] irfsquared
    compute ones=%const(1.0)
    source forcedfactor.src
    *
    compute accept=0
    infobox(action=define,progress,lower=1,upper=n1) 'Monte Carlo Integration'
    do draws=1,n1
       *
       * Make a draw from the posterior for the VAR and compute its impulse
       * responses.
       *
       compute sigmad  =%ranwisharti(svt,wishdof)
       compute p       =%decomp(sigmad)
       compute ranc    =%ran(1.0)
       compute betau   =sxx*ranc*tr(p)
       compute betadraw=betaols+betau
       compute %modelsetcoeffs(varmodel,betadraw)
       *
       * This is changed to unit shocks rather than orthogonalized shocks.
       *
       impulse(noprint,model=varmodel,decomp=%identity(6),results=impulses,steps=nstep)
       gset irfsquared 1 1     = %xt(impulses,t).^2
       gset irfsquared 2 nstep = irfsquared{1}+%xt(impulses,t).^2
       *
       * Do the subdraws over the unit sphere. These give the weights on the
       * orthogonal components.
       *
       do subdraws=1,n2
       ************************************************
       * First Set of Restrictions - Unique Impulse Vector
       ************************************************
          compute v1=%ran(1.0),v1=v1/sqrt(%normsqr(v1))
          compute v=i1=p*v1
          do k=1,KMAX+1
              compute ik=%xt(impulses,k)*v
              if ik(1)>0.or.ik(2)>0.or.ik(3)<0.or.ik(6)>0
                 branch 105
          end do k
          *
          * Meets the first restriction
          * Draw from the orthogonal complement of i1 (last five columns of the
          * factor "f").
          *
          @forcedfactor(force=column) sigmad i1 f
          compute v2=%ran(1.0),v2=v2/sqrt(%normsqr(v2))
          compute v=i2=%xsubmat(f,1,6,2,6)*v2
          *****************************************************
          *  Second Set of Restrictions - Demand Shock
          *****************************************************
          do k=1,KMAX+1
             compute ik=%xt(impulses,k)*v
             if ik(1)<0.or.ik(2)<0.or.ik(3)<0.or.ik(6)<0
                branch 105
          end do k
          *
          * Meets both restrictions
          *
          compute accept=accept+1
          dim goodrespa(accept)(nstep,nvar) goodfevda(accept)(nstep,nvar)
          dim goodresp(accept)(nstep,nvar) goodfevd(accept)(nstep,nvar)
          ewise goodresp(accept)(i,j)=(ik=%xt(impulses,i)*i1),ik(j)
          ewise goodfevd(accept)(i,j)=(ik=(irfsquared(i)*(i1.^2))./(irfsquared(i)*ones)),ik(j)
          ewise goodrespa(accept)(i,j)=(ik=%xt(impulses,i)*i2),ik(j)
          ewise goodfevda(accept)(i,j)=(ik=(irfsquared(i)*(i2.^2))./(irfsquared(i)*ones)),ik(j)
          if accept>=1000
             break
       :105
       end do subdraws
       if accept>=1000
          break
       infobox(current=draws)
    end do draws
    infobox(action=remove)
    *
    * Post-processing. Graph the mean of the responses along with the 16% and 84%-iles
    *
    clear upper lower resp
    *
    spgraph(vfields=3,hfields=2,hlabel='Impulse Responses with Pure-Sign Approach',subhea='Subhea')
    do i=1,nvar
       compute minlower=maxupper=0.0
       smpl 1 accept
       do k=1,nstep
          set work = goodresp(t)(k,i)
          compute frac=%fractiles(work,||.16,.84||)
          compute lower(k)=frac(1)
          compute upper(k)=frac(2)
          compute resp(k)=%avg(work)
          display resp(k)
       end do k
    *
       smpl 1 nstep
       graph(pattern,ticks,number=0,picture='##.##',header='Impulse Responses for '+vl(i)) 3
       # resp
       # upper / 2
       # lower / 2
    end do i
    *
    spgraph(done)
    smpl
    *
    clear upper lower resp
    *
    spgraph(vfields=3,hfields=2,hlabel='SECOND Impulse Responses with Pure-Sign Approach',subhea='Subhea')
    do i=1,nvar
       compute minlower=maxupper=0.0
       smpl 1 accept
       do k=1,nstep
          set work = goodrespa(t)(k,i)
          compute frac=%fractiles(work,||.16,.84||)
          compute lower(k)=frac(1)
          compute upper(k)=frac(2)
          compute resp(k)=%avg(work)
          display resp(k)
       end do k
    *
       smpl 1 nstep
       graph(pattern,ticks,number=0,picture='##.##',header='Impulse Responses for '+vl(i)) 3
       # resp
       # upper / 2
       # lower / 2
    end do i
    *
    spgraph(done)
    smpl
    *
    clear upper lower resp
    *
    spgraph(vfields=3,hfields=2,hlabel='Fraction of Variance Explained with Pure-Sign Approach',subhea='Subhea')
    do i=1,nvar
       compute minlower=maxupper=0.0
       smpl 1 accept
       do k=1,nstep
          set work = goodfevd(t)(k,i)
          compute frac=%fractiles(work,||.16,.50,.84||)
          compute lower(k)=frac(1)
          compute upper(k)=frac(3)
          compute resp(k)=frac(2)
          display resp(k)
       end do k
    *
       smpl 1 nstep
       graph(pattern,ticks,number=0,min=0.0,picture="##.##",header="Fraction Explained for "+vl(i)) 3
       # resp
       # upper / 2
       # lower / 2
    end do i
    *
    spgraph(done)
    smpl
    *
    clear upper lower resp
    *
    spgraph(vfields=3,hfields=2,hlabel='SECOND Fraction of Variance Explained with Pure-Sign Approach',subhea='Subhea')
    do i=1,nvar
       compute minlower=maxupper=0.0
       smpl 1 accept
       do k=1,nstep
          set work = goodfevda(t)(k,i)
          compute frac=%fractiles(work,||.16,.50,.84||)
          compute lower(k)=frac(1)
          compute upper(k)=frac(3)
          compute resp(k)=frac(2)
          display resp(k)
       end do k
    *
       smpl 1 nstep
       graph(pattern,ticks,number=0,min=0.0,picture="##.##",header="Fraction Explained for "+vl(i)) 3
       # resp
       # upper / 2
       # lower / 2
    end do i
    *
    spgraph(done)
    smpl

[TomDoan]

Hi~

Thank you for your instruction.
After removing $, the code seems to wotk.
But I am still confused about how to set up the sign restriction after putting variable1,2 to be equal.
This modification seems only make variable 1,2 to be equal for one period, how to make it last for four periods?
Hope to have your kind responses.
Thanks a lot!

This (the alpha weights on the "swish" matrix) will create responses for the first two variables which are equal to each other for the first four periods. Note that this severely restricts the possible shocks, so it may be difficult for the draws for shocks to satisfy any other constraints.

Code: Select all

   impulse(noprint,model=varmodel,decomp=swish,results=impulses,steps=nstep)
   compute step0=%xt(impulses,1)
   compute step1=%xt(impulses,2)
   compute step2=%xt(impulses,3)
   compute step3=%xt(impulses,4)
   dec rect r(4,6)
   compute %psubmat(r,1,1,||1.0,-1.0,0.0,0.0,0.0,0.0||*step0)
   compute %psubmat(r,2,1,||1.0,-1.0,0.0,0.0,0.0,0.0||*step1)
   compute %psubmat(r,3,1,||1.0,-1.0,0.0,0.0,0.0,0.0||*step2)
   compute %psubmat(r,4,1,||1.0,-1.0,0.0,0.0,0.0,0.0||*step3)
   compute rfperp=%perp(r)
   compute alpha=rfperp*%ranmat(%cols(rfperp),1)
   compute alpha=alpha/sqrt(%normsqr(alpha))