*
* @DLMIRF( options )
* computes and graphs impulse responses and error decompositions for the
* states in a state space model. The transition equation takes the form
* Y(t)=AY(t-1)+Fw(t) where F is an nxp matrix, where p is the number of
* fundamental shocks. The covariance matrix of w is assumed to be the
* identity. To get a different scaling for the w's, post-multiply the F
* for the state space model with a factor of the covariance matrix for
* w(t).
*
* Options:
*   A=NxN matrix A (N=# of states)
*   F=NxP matrix F (p=# of fundamental shocks)
*   STEPS=# of steps [48]
*
*   SHOCKLABELS=VECT[STRINGS] of labels for the shocks ["shock1","shock2",...]
*   VARLABELS=VECT[STRING] of labels for the variables ["var1","var2",...]
*
*   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] computes and displays the error decomposition.
*
*   RESULTS=RECT[SERIES] of responses. The rows are correspond to variables,
*      the columns to shocks.
*
*   [GRAPH]/NOGRAPH. Use NOGRAPH to suppress the graphs.
*
*
* Revision Schedule:
*  05/2007 Written by Tom Doan, Estima.
*  02/2009 Minor change to deal with flat lines
*  03/2009 Add FOOTER option
*  09/2009 Change options to match MCGraphIRF
*  04/2010 Add RESULTS and GRAPH/NOGRAPH options
*  09/2010 Add local report and title handling
*
procedure DLMIRF
option rect          a
option rect          f
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 rect[series]  *results
option switch        graph      1
option string        title
*
local rect[series]  impulses
local vect[series]  resp
local vect[strings] xlabel ylabel
local vect          sumsq impit x0
local rect          impt
local integer 		  i j row step nshocks nshow nstate
local integer       pages nvert nhorz pagetoshow
local integer       shocklower shockupper varlower varupper
local real          norange minlower maxupper
local string        gheader
local vect[integer] grparms
local string        ltitle
local report        dreport
*
if .not.%defined(a).or..not.%defined(f) {
   disp "***DLMIRF needs A and F options to describe model"
   return
}
if %rows(a)<>%cols(a) {
   disp "***DLMIRF - A matrix is" %rows(a) " x " %cols(a) ". Need NxN"
   return
}
if %rows(a)<>%rows(f) {
   disp "**DLMIRF - F matrix is" %rows(f) " x " %cols(f) ". Need " nstate " rows to match A"
   return
}
compute nstate=%rows(a)
compute nshocks=%cols(f)
*
if %defined(varlabels) {
   compute nshow=%size(varlabels)
   compute ylabel=varlabels
}
else {
   compute nshow=nstate
   dim ylabel(nshow)
   ewise ylabel(i)="Variable "+i
}

if %defined(shocklabels)
   compute xlabel=shocklabels
else {
   dim xlabel(nshocks)
   ewise xlabel(i)="Shock "+i
}

dim impulses(nshow,nshocks)
if %defined(results)
   dim results(nshow,nshocks)
do i=1,nshocks
   do j=1,nshow
      set impulses(j,i) 1 steps = 0.0
   end do j
   compute x0=%xcol(f,i)
   do step=1,steps
      do j=1,nshow
         compute impulses(j,i)(step)=x0(j)
      end do j
      compute x0=a*x0
   end do step
   if %defined(results) {
      do j=1,nshow
         set results(j,i) 1 steps = impulses(j,i)
      end do j
   }
end do i
*
if graph==0
   goto skipgraph

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(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
            compute gheader=""
			if page==1 {
            graph(ticks,min=minlower,max=maxupper,number=0,header=gheader) 1 j i
            # resp(j)  1 steps
         }
         else {
            graph(ticks,min=minlower,max=maxupper,number=0,header=gheader) 1
            # resp(j)  1 steps
			}
      end do j
   end do i
   if page<>2
      spgraph(done)
end do pagetoshow
grparm(recall=grparms)

:skipgraph

if errors {
   if %defined(title)
      compute ltitle=title
   else
      compute ltitle="State-Space Decomposition of Variance"

   report(use=dreport,action=define,title=ltitle)
   compute row=1
   do i=1,nshow
      report(use=dreport,atrow=row,atcol=1,span) "Decomposition of Variance, Series "+ylabel(i)
      compute row=row+1
		report(use=dreport,atrow=row,atcol=1) "Step" xlabel
      compute sumsq=%zeros(nshocks,1)
      do step=1,steps
         compute impt =%xt(impulses,step)
         compute impit=%xrow(impt,i)
         compute sumsq=sumsq+impit.*impit
         report(use=dreport,atrow=row,atcol=1,fillby=rows) step sumsq/%sum(sumsq)*100.0
         compute row=row+1
      end do step
      compute row=row+1
   end do i
   report(use=dreport,action=format,picture="*.###")
   report(use=dreport,action=show,window=ltitle)
}

end DLMIRF

