The RATS Software Forum

Posted: **Tue Aug 01, 2017 8:30 am**

VWEIGHTS would just be V1 (not V1~other shock descriptors)

Posted: **Tue Aug 01, 2017 9:01 am**

Here is the adjusted code I wrote to save the shocks:

open data uhligdata.xls
calendar(m) 1965:1
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
*


source uhligfuncs.src



*** VAR Setup and Storage Matrices

compute lags=12                     ;*Number of lags
compute nstep=60                    ;*Number of response steps
compute n1= 10000                   ;*Number of draws for Gibbs-Loop
compute n2= 50000                   ;*Number of draws from the unit sphere for each VAR draw
compute nburn = 9000                ;*Number of burn-in-draws
compute keep = n1 - nburn           ;*Number of kept draws
compute nvari = 6                   ;*Number of variables
compute kmin = 1                    ;*Number of Steps constrained (Lower Bouund)
compute kmax = 5                    ;*Number of Steps constrained (Upper Bouund)
compute ncoef = nvari*lags          ;*Number of coefficients per equation
compute nkeep = 10000               ;*Number of excepted draws that we actually want


compute tight=.1
compute other=.5


system(model=varmodel)
variables gdpc1 gdpdef cprindex fedfunds bognonbr totresns
lags 1 to 12
end(system)

estimate(noprint,ols,sigma)                                       ;*initial set of estimations
impulse(noprint,model=varmodel,steps=nstep,results=impulses)      ;*initial set of of impulse responses

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


compute nvar=%nvar                                                ;* Number of equations estimated
dec vect olssee(nvar)                                             ;* vector of variances
ewise olssee(i)=sqrt(%sigma(i,i))                                 ;* standard deviations for rescaling

declare series[rect] bgibbs                                       ;* Storage Series of Rectangulars to save the VAR Coefficient Draws

gset bgibbs 1 keep = %zeros(ncoef,nvari)                          ;* To dimension the storage matrix inside the series

compute hstart=%regstart(),hend=%regend()                         ;* Save the estimation range for computing residuals

declare vect[series] shock(nkeep)

***********************************
********* Minnesota Prior *********
***********************************

source BVARBuildPriorMN.src

@BVARBuildPriorMN(model=varmodel,tight=tight,other=other) HBPRIORS HPRIORS
@BVARFinishPrior hbpriors hpriors bprior hprior

compute sigmad=%sigma
cmom(noprint,model=varmodel)



*********************************
********* Gibbs Sampler *********
*********************************

declare vect[rect] goodresp(nkeep) goodfevd(nkeep)
declare vector ik a(nvar) ones(nvar)
declare series[rect] irfsquared
compute ones=%const(1.0)                                          ;* Fills matrix with a constant value

compute accept = 0

infobox(action=define,progress,lower=1,upper=n1) "Gibbs Sampler"
do draw= 1,n1
   infobox(current=draw)
   *
   * Draw b given sigma
   *
   compute bdraw =%ranmvkroncmom(%cmom,inv(sigmad),hprior,bprior)
   compute rssmat=%sigmacmom(%cmom,bdraw)
   *
   * Draw sigma given b
   *
   compute sigmad = %ranwisharti(%decomp(inv(rssmat)),%nobs)
   compute swish = %decomp(sigmad)
   if draw<=nburn
      next

   compute %modelsetcoeffs(varmodel,bdraw)
   forecast(model=varmodel,from=hstart,to=hend,errors=resids,static)
   impulse(noprint,model=varmodel,results=impulses,steps=nstep,decomp=swish)

   gset irfsquared 1 1 = %xt(impulses,t).^2
   gset irfsquared 2 nstep = irfsquared{1}+%xt(impulses,t).^2


      do subdraws=1,n2
         compute a=%ransphere(nvar)
         *
         * 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).
         *
         if UhligAccept(a,KMIN,KMAX+1,||+4,-3,-2,-5||)==0
         goto reject
         *
         * This is an accepted draw. Copy the information out. If we
         * have enough good ones, drop out of the overall loop.
         *
         compute accept=accept+1

         dim goodresp(accept)(nstep,nvari)
         dim goodfevd(accept)(nstep,nvari)

         ewise goodresp(accept)(i,j)=ik=%xt(impulses,i)*a,ik(j)                              ;*%xt(impulses,k) pulls the [nvar x nvar] matrix of responses to the Cholesky shocks at step k.
         ewise goodfevd(accept)(i,j)=ik=(irfsquared(i)*(a.^2))./(irfsquared(i)*ones),ik(j)

         if accept>=nkeep
         break
         :reject
      end do subdraws

   dim shock(accept)
   set shock(accept) = %dot(a,inv(sigmad)*%xt(resids,t))

   if accept>=nkeep
   break

   * Bookkeeping of VAR coefficients:

   compute bgibbs(draw-nburn) = bdraw                                                       ;* Saves the Draws of the VAR Coefficients

end do draw
infobox(action=remove)


graph 1
#shock(2)

when I want to graph one of the saved shocks I get the following error message: "## SR14. Empty Range on Series SHOCK(2)."

Can you help?

Thanks

Best

Jules

Posted: **Tue Aug 01, 2017 9:41 am**

You made major (and incorrect) changes to this rather than just reducing nshocks to 1 and fixing VWEIGHTS. Go back to the original and keep the basic structure.

dim shocks(accept)(nshocks)
      compute vweights=v1~v2~v3r~v3e
      do j=1,nshocks
         set shocks(accept)(j) = %dot(%xcol(vweights,j),inv(sigmap)*%xt(resids,t))
      end do j
      if accept>=nkeep
         break

Posted: **Wed Aug 02, 2017 10:44 am**

You are right! Here is a corrected code:

Code: Select all

open data uhligdata.xls
calendar(m) 1965:1
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
*


source uhligfuncs.src



*** VAR Setup and Storage Matrices

compute lags=12                     ;*Number of lags
compute nstep=60                    ;*Number of response steps
compute n1= 10000                   ;*Number of draws for Gibbs-Loop
compute n2= 50000                   ;*Number of draws from the unit sphere for each VAR draw
compute nburn = 9000                ;*Number of burn-in-draws
compute keep = n1 - nburn           ;*Number of kept draws
compute nvari = 6                   ;*Number of variables
compute kmin = 1                    ;*Number of Steps constrained (Lower Bouund)
compute kmax = 5                    ;*Number of Steps constrained (Upper Bouund)
compute ncoef = nvari*lags          ;*Number of coefficients per equation
compute nkeep = 10000               ;*Number of excepted draws that we actually want


compute tight=.1
compute other=.5


system(model=varmodel)
variables gdpc1 gdpdef cprindex fedfunds bognonbr totresns
lags 1 to 12
end(system)

estimate(noprint,ols,sigma)                                       ;*initial set of estimations
impulse(noprint,model=varmodel,steps=nstep,results=impulses)      ;*initial set of of impulse responses

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


compute nvar=%nvar                                                ;* Number of equations estimated
dec vect olssee(nvar)                                             ;* vector of variances
ewise olssee(i)=sqrt(%sigma(i,i))                                 ;* standard deviations for rescaling

declare series[rect] bgibbs                                       ;* Storage Series of Rectangulars to save the VAR Coefficient Draws

gset bgibbs 1 keep = %zeros(ncoef,nvari)                          ;* To dimension the storage matrix inside the series

compute hstart=%regstart(),hend=%regend()                         ;* Save the estimation range for computing residuals

declare vect[vect[series]] shocks(nkeep)

***********************************
********* Minnesota Prior *********
***********************************

source BVARBuildPriorMN.src

@BVARBuildPriorMN(model=varmodel,tight=tight,other=other) HBPRIORS HPRIORS
@BVARFinishPrior hbpriors hpriors bprior hprior

compute sigmad=%sigma
cmom(noprint,model=varmodel)



*********************************
********* Gibbs Sampler *********
*********************************

declare vect[rect] goodresp(nkeep) goodfevd(nkeep)
declare vector ik a(nvar) ones(nvar)
declare series[rect] irfsquared
compute ones=%const(1.0)                                          ;* Fills matrix with a constant value

compute nshocks=1                                                 ;* Number of shocks
dec vect[vect[series]] shocks(nkeep)                              ;* Storage vector for the series of shocks


compute accept = 0

infobox(action=define,progress,lower=1,upper=n1) "Gibbs Sampler"
do draw= 1,n1
   infobox(current=draw)
   *
   * Draw b given sigma
   *
   compute bdraw =%ranmvkroncmom(%cmom,inv(sigmad),hprior,bprior)
   compute rssmat=%sigmacmom(%cmom,bdraw)
   *
   * Draw sigma given b
   *
   compute sigmad = %ranwisharti(%decomp(inv(rssmat)),%nobs)
   compute swish = %decomp(sigmad)
   if draw<=nburn
      next

   compute %modelsetcoeffs(varmodel,bdraw)
   forecast(model=varmodel,from=hstart,to=hend,errors=resids,static)
   impulse(noprint,model=varmodel,results=impulses,steps=nstep,decomp=swish)

   gset irfsquared 1 1 = %xt(impulses,t).^2
   gset irfsquared 2 nstep = irfsquared{1}+%xt(impulses,t).^2


      do subdraws=1,n2
         compute a=%ransphere(nvar)
         *
         * 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).
         *
         if UhligAccept(a,KMIN,KMAX+1,||+4,-3,-2,-5||)==0
         goto reject
         *
         * This is an accepted draw. Copy the information out. If we
         * have enough good ones, drop out of the overall loop.
         *
         compute accept=accept+1

         dim goodresp(accept)(nstep,nvari)
         dim goodfevd(accept)(nstep,nvari)

         ewise goodresp(accept)(i,j)=ik=%xt(impulses,i)*a,ik(j)                              ;*%xt(impulses,k) pulls the [nvar x nvar] matrix of responses to the Cholesky shocks at step k.
         ewise goodfevd(accept)(i,j)=ik=(irfsquared(i)*(a.^2))./(irfsquared(i)*ones),ik(j)

         dim shocks(accept)(nshocks)
         compute vweights=a
         do j=1,nshocks
            set shocks(accept)(j) = %dot(%xcol(vweights,j),inv(swish)*%xt(resids,t))
         end do j
         if accept>=nkeep
            break

         :reject
      end do subdraws



   if accept>=nkeep
   break

   * Bookkeeping of VAR coefficients:

   compute bgibbs(draw-nburn) = bdraw                                                       ;* Saves the Draws of the VAR Coefficients

end do draw
infobox(action=remove)



graph 3
# shocks(2)(1)
# shocks(1)(1)
# shocks(10000)(1)

That should be correct right?

Just one more question. Is it possible to delete a series after it is loaded/generated?
Something like

delte shocks(10000)(1)

or so?

Because I saved the 100000 shock series and after I did some analysis to them, I need some space for further analysis, thats why I want to delte them.

Best Jules

Posted: **Wed Aug 02, 2017 5:10 pm**

dim shocks(0)

should do it.

The RATS Software Forum

BVAR Minnesota Prior IRF

Re: BVAR Minnesota Prior IRF

Re: BVAR Minnesota Prior IRF

Re: BVAR Minnesota Prior IRF

Re: BVAR Minnesota Prior IRF

Re: BVAR Minnesota Prior IRF