*
* @VARIRF(model=VAR model,other options)
*
* VARIRF organizes the graphs of an impulse response function from an already
* estimated VAR.
*
* Options:
*   MODEL=VAR model whose IRF's are to be calculated
*   FACTOR=Factor of the covariance matrix [Cholesky of %SIGMA]
*   STEPS=number of forecast steps[48]
*
*   SHOCKLABELS=VECT[STRINGS] of labels for the shocks [dependent variables]
*   VARLABELS=VECT[STRING] of labels for the variables [dependent variables]
*
*   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.
*
*   ERRORS/[NOERRORS] Produce a forecast error variance decomposition
*   ACCUMULATE=||vector of positions||. Indicates which variables need to be
*    accumulated (integrated)
*
*   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.
*
*   PICTURE="picture for value labels" [automatically chosen for each series]
*
*   NUMBER=labeling number assigned to first entry plotted [0]
*
* Revision Schedule:
*   07/2005 Written by Tom Doan, Estima.
*   05/2009 Add FACTOR option for compatibility with new syntax
*   05/2009 Add ACCUMULATE option
*   10/2009 Change options to match similar procedures
*   10/2017 Add INCLUDE option
*   08/2024 Add PICTURE option
*   06/2025 Add NUMBER option
*   06/2025 Adapt to allow for a non-square FACTOR
*
procedure VARIRF
*
option model         model
option integer       steps      48
option string        header
option string        footer
option choice        page       1   all  one  byshock   byvariable
option integer       columns
option switch        errors     0
option vect[strings] shocklabels
option vect[strings] varlabels
option vect[integer] accumulate
option vect[integer] include
option switch        print      0
option string        picture
option integer       number     0
*
* These are synonyms
*
option rect          factor
option rect          decomp
*
local rect[series]   impulses
local vect[series]   resp
local vect[strings]  ylabel
local vect[strings]  xlabel
local vect[integer]  depvars yshow
local integer        neqn
local integer        i j
local string         header
local rect           factord
local integer        pages nvert nhorz
local integer        ntargets nshocks nshow
local vect[integer]  grparms
local integer        pagetoshow shocklower shockupper varlower varupper
local real           minlower maxupper norange
local string         gheader
*
if .not.%defined(model) {
   display "Syntax: @VARIRF(model=VAR model,steps=# of steps,...)"
   return
}
compute neqn=%modelsize(model)
compute ntargets=neqn
compute depvars=%modeldepvars(model)
*
* Use the varlabels if provided, else use the labels of the dependent
* variables.
*
if %defined(include)
   compute nshow=%size(include),yshow=include
else
   compute nshow=ntargets,yshow=%seq(1,nshow)

dim ylabel(nshow)
if %defined(varlabels)
   ewise ylabel(i)=varlabels(yshow(i))
else
   ewise ylabel(i)=%l(depvars(yshow(i)))
*
* Use the FACTOR or DECOMP provided, otherwise take the standard
* factorization of the %SIGMA matrix.
*
if %defined(factor)
   compute factord=factor
else
if %defined(decomp)
   compute factord=decomp
else
   compute factord=%decomp(%sigma)
*
compute nshocks=%cols(factord)
*
* Use the shocklabels if provided, otherwise label the shocks the same
* as the responses.
*
if %defined(shocklabels)
   compute nshocks=%imin(nshocks,%size(shocklabels))

dim xlabel(nshocks)

if %defined(shocklabels)
   ewise xlabel(i)=shocklabels(i)
else
   ewise xlabel(i)=%l(depvars(i))
*
* This computes the full set of impulse responses, which are in the
* series in impulses. impulses(i,j) is the response of variable i to a
* shock in j.
*
dim impulses(neqn,nshocks)
impulse(model=model,result=impulses,print=print,decomp=factord,steps=steps)
*
* Accumulate the requested responses
*
if %defined(accumulate) {
   dofor i = accumulate
      if i<=0.or.i>neqn
         next
      do j=1,nshocks
         acc impulses(i,j) 1 steps
      end do j
   enddofor i
}
*
* This loop graphs the responses of a variable to all shocks.
* These don’t have to be normalized.
*
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
   }

*
* Save the graph parms and reset the sizes
*

compute grparms = %grparm()
grparm(bold) hlabel 18 matrixlabel 14
grparm  axislabel 24

dim resp(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 = impulses(yshow(i),j)
         compute maxupper=%max(maxupper,%maxvalue(resp(j)))
         compute minlower=%min(minlower,%minvalue(resp(j)))
      end do j
      *
      * Adjust if responses seem to be roundoff error
      *
      compute norange=%noprec(.01*(maxupper-minlower))
      if norange.and.%noprec(.01*minlower)
         compute maxupper=1.0
      else
      if norange.and.maxupper>0.0
         compute minlower=0.0
      else
      if norange.and.minlower<0.0
         compute maxupper=0.0
      do j=shocklower,shockupper
         if page==3
            compute gheader=ylabel(i)
         else
         if page==4
            compute gheader=xlabel(j)
         else
         if page==2
            compute gheader="Responses of "+ylabel(i)+" to "+xlabel(j)
         else
            compute gheader=""
         if page==1 {
            graph(ticks,min=minlower,max=maxupper,header=gheader,picture=picture,number=number,col=j,row=i)
            # resp(j)  1 steps
         }
         else {
            graph(ticks,min=minlower,max=maxupper,header=gheader,picture=picture,number=number) 1
            # resp(j)  1 steps
         }
      end do j
   end do i
   if page<>2
      spgraph(done)
end do pagetoshow
grparm(recall=grparms)
*
if errors
   errors(labels=shocklabels,model=model,steps=steps,$
     window="Decomposition of Variance",decomp=factor)
end varirf
