The RATS Software Forum

The attached zip file has replication files for Mountford & Uhlig(2009), "What are the Effects of Fiscal Policy Shocks?", Journal of Applied Econometrics (to appear). mu2009b1.prg does the analysis with 4 period sign constrained shocks, including the combination shocks to do the various forms of fiscal policy experiments. mu2009b2.prg does the analysis with "delayed" responses. There's a common source file with the set up code for the model.

To (greatly) simplify the code, we wrote a function which computes the "penalty function" for responses based upon the sign constraints.

Code: Select all

*
* UhligPenalty(q,first,last,constrained) returns the penalty function
* for the weight vector <<q>> (applied to the impulse responses in the
* global VECT[SERIES] impulses). <<first>> and <<last>> are the range of
* responses constrained (1 = impact response). <<constrained>> is a
* VECT[INTEGER] with the variable positions being constrained. Use +slot
* for a positivity constraint and -slot for a negativity constraint.
* That is, if you want the first 4 responses to be positive on the 2nd
* variable and negative on the third, you would use the function
*
* UhligPenalty(q,1,4,||2,-3||)
*
function UhligPenalty q first last constrained
type real      UhligPenalty
type vector    q
type integer   first last
type vect[int] constrained
*
local integer i k
local real    func
local vector  ik
*
compute func=0.0
do k=first,last
   compute ik=(%xt(impulses,k)*q)./scales
   do i=1,%rows(constrained)
      if constrained(i)<0
         compute value= ik(-constrained(i))
      else
         compute value=-ik(constrained(i))
      compute func=func+%if(value<0.0,value,100*value)
   end do i
end do k
compute UhligPenalty=func
end

This also uses a different method for handling the minimization of the penalty function than was done in the program for the Uhlig JME 2005 paper. Instead of doing each optimization twice, then rejecting draws where the modes don't match up, it does a slower broad scan with the genetic algorithm as the preliminary method, which seems in practice to find the global optimum with very high reliability. With this (fairly large) model, it's able to do 250 draws in about five minutes, which is much faster than the Gauss code provided by the authors.

Dear Tom,

I can't find hte mcgraphirf.src in the zipped file.
Would you please tell me how to get it?
thank you very much!

It's on the web site in the procedures page. It's also on the forum at

viewtopic.php?f=7&t=332

Note that there's a "search" button in the top right which will search only the forum.

Dear Tom,

Thank you very much.
But mu2009b2.prg needs to add one line: source mujaesetup.src.

Please please instruct me
1.how to modify mu2009b1.prg to compute forecast error variance decomposition?
2.how to modify mcgraphirf.src to show the sign restrictions on the impulse responses?

Best regards.

I'm not sure what you mean by 2.

For example, the government spending shock requires the government spending increases for 4 periods, then you will see a line on the period 4 in
impulse response of government spending.

That would be done with a GRID=(T==KMAX) or a SHADING=(T<=KMAX) option on GRAPH.

The IRFSQUARED's are already being computed in the code, so you would just add the same code as is done in the Uhlig JME replication to do the variance decomposition.

Hi Tom,

Thanks for your reply.
But I still can't figure out how to add GRID=(T==KMAX) or a SHADING=(T<=KMAX) on the graph into the code_mcgraphirf.src as followed.
especially when we just wanna grid the response of variable 2 in the impulse response graph of government spending shock, grid the response of variable 1, 3, 8, 9 in the impulse response graph of business cycle shock. how can we make it?

Code: Select all

*
* @MCGraphIRF(model=model used in generating responses,other options )
*
* Graphs error bands for impulse response functions using the
* information already computed in %%responses. These can be calculated
* using @MCVARDoDraws or some similar procedure.
*
* %%responses should have the following structure:
*   the number of draws is NDRAWS,
*   the number of steps is NSTEPS
*   the number of variables is NVAR and
*   the number of shocks is NSHOCKS.
*   (NSHOCKS often is equal to NVAR, but doesn't have to be).
*
* %%responses is a VECT[RECT] with outer dimensions NDRAWS. Each draw is
* represented by a RECT array, with NVAR*NSHOCKS rows and NSTEPS
* columns. The rows are blocked by shocks so the first set of NVAR
* elements in a column are the responses to the first shock, the second
* set are the responses to the second shock. If you %VEC a set of
* impulse responses produced by IMPULSES, that's how they will be
* blocked.
*
* Options:
*   MODEL=model used in generating responses
*   SHOCKLABELS=VECT[STRINGS] of labels for the shocks [dependent variables of model]
*   VARLABELS=VECT[STRING] of labels for the variables [dependent variables of model]
*
*   INCLUDE=||list of dependent variables to show (in order, by position in model)|| [all]
*      INCLUDE can be used to eliminate some variables which aren't of
*      special interest or to rearrange the panes on a graph. For instance,
*      in a 6 variable VAR, INCLUDE=||6,1,3,4|| will show only the responses
*      of those four variables (based upon their positions in the model) and
*      will put the 6th variable in the top left.
*
*   CENTER=[MEAN]/MEDIAN/INPUT
*   IMPULSES=RECT[SERIES] of central impulse responses for the graph
*      This chooses what is shown as the "estimate" of the IRF. MEAN and
*      MEDIAN are the mean and median of the drawn responses. INPUT means
*      that you are including an IMPULSES option with RECT[SERIES] which are
*      to be used. This needs to be in the format created by the IMPULSE
*      instruction.
*
*   PERCENTILES=||percentiles for lower and upper bounds|| [||.16,.84||]
*   STDDEV=# of standard deviations from mean for lower and upper bounds [not used]
*      STDDEV is used for doing error bands based upon multiples of the
*      sample standard deviations. For instance, STDDEV=1.0 will give upper
*      and lower bounds that are one standard deviation above and below the
*      central response. PERCENTILES is the default.
*
*   HEADER=' title for graph '
*   FOOTER=' footer for graph '
*   PAGE=[ALL]/ONE/BYSHOCK/BYVARIABLE
*   COLUMNS=# of columns on a page [depends upon number of graphs]
*      The PAGE option selects the layout of a single page of graphs.
*
*      PAGE=ALL (the default) does all responses to all shocks on a single
*        page (shocks in columns, dependent variables in rows). Note that the
*        panes get a bit small when you get more than four variables.
*
*      PAGE=ONE does one combination of shock and variable per page.
*
*        PAGE=BYSHOCK does a separate page for each shock, with all responses
*        arranged in one or more columns.
*
*        PAGE=BYVARIABLE does a separate page for each variable, with all
*        shocks arranged in one or more columns.
*
* Revision Schedule:
*   06/2009 Stripped from montevar.src procedure. Options added.
*
procedure MCGraphIRF
*
option model         model
option vect[strings] shocklabels
option vect[strings] varlabels
option vector        percentiles  ||.16,.84||
option real          stddev
option choice        center       1   mean median input
option rect[series]  impulses
option string        header
option string        footer
option choice        page         1   all  one  byshock   byvariable
option integer       columns
option vect[integer] include
*
local  vect[integer] grparms
local  vect[strings] xlabel ylabel
local  vect[integer] depvars yshow
local  integer       draws steps nshocks nshow npercent
local  vect[series]  upper lower upperx lowerx resp
local  real          minlower maxupper sigma
local  integer       i j k
local  vector        work
local  vector        frac
local  integer       ntargets pages nvert nhorz pagetoshow
local  integer       shocklower shockupper varlower varupper
local  string        gheader
local  vect          request
*
declare vect[rect]   %%responses
*
if %rows(%%responses)==0 {
   disp "####@MCGraphIRF - needs procedure to draw responses first (%%responses not defined)"
   return
}

if .not.%defined(model) {
   disp "Syntax: @MCGraphIRF(model=model used in generating responses,other options)"
   return
}

if center==3.and..not.%defined(impulses) {
   disp "####@MCGraphIRF(CENTER=INPUT) needs IMPULSES option"
   return
}
*
* Take the number of variables (targets) out of the model
*
compute ntargets=%modelsize(model)
*
* Hack the dimensions of everything else out of the <<%%responses>> array
*
compute draws=%rows(%%responses)
compute steps=%cols(%%responses(1))
compute nshocks=%rows(%%responses(1))/ntargets
*
dim work(draws)
*
* Set up the labels for the matrix of graphs
*
compute depvars=%modeldepvars(model)

if %defined(include)
   compute nshow=%rows(include),yshow=include
else
   compute nshow=ntargets,yshow=%seq(1,nshow)

dim xlabel(nshocks) ylabel(nshow)

if %defined(shocklabels)
   ewise xlabel(i)=shocklabels(i)
else
   ewise xlabel(i)=%l(depvars(i))

if %defined(varlabels)
   ewise ylabel(i)=varlabels(yshow(i))
else
   ewise ylabel(i)=%l(depvars(yshow(i)))

*
* Save the graph parms and reset the sizes
*

compute grparms = %grparm()
grparm(bold) hlabel 18 matrixlabels 14
grparm  axislabel 24
*
* Figure out how many percentiles we need, combining the requests from
* the bands (in percentiles option) and center (if median).
*
dim request(0)
compute npercent=0
if %defined(percentiles).and..not.%defined(stddev)
   compute npercent=%rows(percentiles)

if center==2 {
   if npercent>0
      compute request=percentiles~~||.50||
   else
      compute request=||.50||
}
else
if npercent>0
   compute request=percentiles

if page==1
   compute pages=1,nvert=nshow,nhorz=nshocks
else
if page==2
   compute pages=nshow*nshocks,nvert=nhorz=1
else
if page==3 {
   compute pages=nshocks
   if %defined(columns).and.columns>0
      compute nhorz=columns
   else
      compute nhorz=fix(sqrt(nshow))
   compute nvert=(nshow-1)/nhorz+1
   }
else
if page==4 {
   compute pages=nshow
   if %defined(columns).and.columns>0
      compute nhorz=columns
   else
      compute nhorz=fix(sqrt(nshocks))
   compute nvert=(nshocks-1)/nhorz+1
   }

dim upper(nshocks) lower(nshocks) resp(nshocks)
if npercent>2
   dim upperx(nshocks) lowerx(nshocks)

do pagetoshow=1,pages
   if page==1 {
      spgraph(header=header,footer=footer,xpos=both,xlab=xlabel, $
        ylab=ylabel,vlab="Responses of",vfields=nshow,hfields=nshocks)
      compute shocklower=1,shockupper=nshocks
      compute varlower=1,varupper=nshow
   }
   else
   if page==3 {
      spgraph(footer="Responses to "+xlabel(pagetoshow),vfields=nvert,hfields=nhorz)
      compute shocklower=pagetoshow,shockupper=pagetoshow
      compute varlower=1,varupper=nshow
   }
   else
   if page==4 {
      spgraph(footer="Responses of "+ylabel(pagetoshow),vfields=nvert,hfields=nhorz)
      compute shocklower=1,shockupper=nshocks
      compute varlower=pagetoshow,varupper=pagetoshow
   }
   else
   if page==2 {
      compute shocklower=(pagetoshow-1)/nshow+1,shockupper=shocklower
      compute varlower=%clock(pagetoshow,nshow),varupper=varlower
   }
   do i=varlower,varupper
      compute minlower=maxupper=0.0
      do j=shocklower,shockupper
         set resp(j)  1 steps = 0.0
         set lower(j) 1 steps = 0.0
         set upper(j) 1 steps = 0.0
         if npercent>2 {
            set lowerx(j) 1 steps = 0.0
            set upperx(j) 1 steps = 0.0
         }
         do k=1,steps
            *
            * Extract the record (across draws) for the response being graphed
            *
            ewise work(t)=%%responses(t)((j-1)*ntargets+yshow(i),k)
            *
            * Compute whatever percentiles are needed
            *
            compute frac=%fractiles(work,request)
            *
            * Choose the central value
            *
            if center==1
               compute resp(j)(k)=%avg(work)
            else
            if center==2
               compute resp(j)(k)=frac(npercent+1)
            else
               compute resp(j)(k)=impulses(yshow(i),j)(k)

            if %defined(stddev) {
               compute sigma=sqrt(%normsqr(work)/draws-resp(j)(k)^2)
               compute lower(j)(k)=resp(j)(k)-stddev*sigma
               compute upper(j)(k)=resp(j)(k)+stddev*sigma
            }
            else {
               if npercent>2 {
                  compute lower(j)(k)=frac(1)
                  compute lowerx(j)(k)=frac(2)
                  compute upperx(j)(k)=frac(3)
                  compute upper(j)(k)=frac(4)
               }
               else {
                  compute lower(j)(k)=frac(1)
                  compute upper(j)(k)=frac(2)
               }
            }
         end do k
         compute maxupper=%max(maxupper,%maxvalue(upper(j)))
         compute minlower=%min(minlower,%minvalue(lower(j)))
      end do j
      *
      do j=shocklower,shockupper
         if page==3
            compute gheader=ylabel(i)
         else
         if page==4
            compute gheader=xlabel(j)
         else
            compute gheader=""
         if page==1.and.npercent>2 {
            graph(ticks,min=minlower,max=maxupper,number=0,header=gheader) 5 j i
            # resp(j)  1 steps
            # upper(j) 1 steps 2
            # lower(j) 1 steps 2
            # upperx(j) 1 steps 3
            # lowerx(j) 1 steps 3
         }
         else
         if page==1 {
            graph(ticks,min=minlower,max=maxupper,number=0,header=gheader) 3 j i
            # resp(j)  1 steps
            # upper(j) 1 steps 2
            # lower(j) 1 steps 2
         }
         else
         if page<>1.and.npercent>2 {
            graph(ticks,min=minlower,max=maxupper,number=0,header=gheader) 5
            # resp(j)  1 steps
            # upper(j) 1 steps 2
            # lower(j) 1 steps 2
            # upperx(j) 1 steps 3
            # lowerx(j) 1 steps 3
         }
         else {
            graph(ticks,min=minlower,max=maxupper,number=0,header=gheader) 3
            # resp(j)  1 steps
            # upper(j) 1 steps 2
            # lower(j) 1 steps 2
         }
      end do j
   end do i
   if page<>4
      spgraph(done)
end do pagetoshow
*
grparm(recall=grparms)
end MCGraphIRF


[b][/b]

The "i" in this loop is the variable:

Code: Select all

do i=varlower,varupper
      compute minlower=maxupper=0.0
      do j=shocklower,shockupper
      ....

If you want a GRID each time series 2 is graphed, add something like:

Code: Select all

do i=varlower,varupper
   if i==2
      set grid 1 steps = t==KMAX
   else
      set grid 1 steps = 0.0


then add the option GRID=GRID to all the GRAPH instructions inside the loop.

Thanks for your code.

However, I have a question that why simplex algorithm is needed to be performed after genetic algorithm? I am curious that why the genetic algorithm cannot be used alone for obtaining the minimum point since genetic algorithm itself is the optimization tool like simplex algorithm.

dennis0125hk wrote:Thanks for your code.

However, I have a question that why simplex algorithm is needed to be performed after genetic algorithm? I am curious that why the genetic algorithm cannot be used alone for obtaining the minimum point since genetic algorithm itself is the optimization tool like simplex algorithm.

Genetic is good for a broad scan of the space, but once you've picked a "hill", it takes much longer to reach final convergence. Since this optimization is inside a loop, efficiency matters.

lhlee0506 wrote:Dear Tom,

Thank you very much.
But mu2009b2.prg needs to add one line: source mujaesetup.src.

Please please instruct me
1.how to modify mu2009b1.prg to compute forecast error variance decomposition?
2.how to modify mcgraphirf.src to show the sign restrictions on the impulse responses?

Best regards.

Dear sir ,
could you give me some suguessitions about how to modify mu2009b1.prg to compute forecast error variance decomposition??
Thank you very much!

dear sir:
In this Mountford$Uhlig JAE 2009 replication files,there are four structural shocks--a business cycle shock ,a monetary policy shock ,a government revenue shock and a government spending shock.
The question is :how to get the four structural shocks' innovations,then I could see the shocks' track directly?
thanks a lot.