*
* @MCFEVDTable(model=model used in generating responses,other options)
*
* Computes error bands for the forecast error variance decomposition for
* a VAR using using the information already computed in %%responses.
* These can be calculated using @MCVARDoDraws or some similar procedure.
*
* Note: this assumes that the shocks used in generating %%responses are
* orthogonal and produce a complete factorization of the covariance
* matrix. It cannot be used with isolated shocks (from, for instance,
* sign restrictions).
*
* %%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.
*
*   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
*     mean response. PERCENTILES is the default.
*
* Revision Schedule:
*   06/2009 Written by Tom Doan, Estima.
*   07/2009 Revised to use different blocking on %%responses
*
procedure MCFEVDTable
option model         model
option vect[integer] horizons     ||1,2,3,4,8,12,40||
option vect[strings] varlabels
option vect[strings] shocklabels
option vector        percentiles  ||.16,.50,.84||
option vect[integer] include
option integer       stddev
option choice        table        1   byvariable  byshock
*
local vect[rect]     fevd
local vect[strings]  xlabel ylabel
local rect           squares
local vect           running totals
local integer        draws ntargets nshocks nshow steps
local integer        i k draw baserow inner outer
local vect           work frac
local vect[integer]  depvars yshow
*
declare vect[rect]   %%responses
*
if %rows(%%responses)==0 {
   disp "####@MCFEVDTable - needs procedure to draw responses first (%%responses not defined)"
   return
}

if .not.%defined(model) {
   disp "Syntax: @MCFEVDTable(model=model used in generating responses,other options)"
   return
}
compute ntargets=%modelsize(model)
*
* Hack the dimensions of everything out of the <<%%responses>> array
*
compute draws=%rows(%%responses)
compute steps=%cols(%%responses(1))
compute nshocks=%rows(%%responses(1))/ntargets
*
* Set up the labels for the tables
*
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)))

local vect[rect] fevd(draws)
*
dim squares(ntargets*nshocks,steps) running(ntargets*nshocks) totals(ntargets)
*
do draw=1,draws
   dim fevd(draw)(ntargets*nshocks,steps)
   *
   * Take elementwise squares of the drawn orthogonalized responses
   *
   compute squares=%%responses(draw).^2
   *
   * Compute the running sums of the squared responses
   *
   compute running=%zeros(ntargets*nshocks,1)
   do k=1,steps
      ewise running(i)=running(i)+squares(i,k)
      *
      * This computes the sums across rows of the squared responses,
      * which will be the total variance for each variable. Dividing the
      * running sums by the total gives the fraction explained by each
      * shock at step <<k>>
      *
      compute totals=%vectorect(running,ntargets)*%fill(nshocks,1,1.0)
      do i=1,ntargets*nshocks
         compute fevd(draw)(i,k)=running(i)/totals(%clock(i,ntargets))
      end do i
   end do k
end do draw
*
* Convert information to table
*
dim work(draws)

if table==1
   compute outer=nshow,inner=nshocks
else
   compute outer=nshocks,inner=nshow

do i=1,outer
   report(action=define)
   if table==1 {
      report(atrow=1,atcol=1,span) "Decomposition of Variance for "+ylabel(i)
      report(atrow=2,atcol=1) "Horizon" xlabel
   }
   else {
      report(atrow=1,atcol=1,span) "Percentage of Variance due to "+xlabel(i)
      report(atrow=2,atcol=1) "Horizon" ylabel
   }
   compute baserow=3
   do k=1,%rows(horizons)
      report(atrow=baserow,atcol=1) horizons(k)
      do j=1,inner
         if table==1
            ewise work(t)=fevd(t)((j-1)*ntargets+yshow(i),horizons(k))
         else
            ewise work(t)=fevd(t)((i-1)*ntargets+yshow(j),horizons(k))
         compute frac=%fractiles(work,percentiles)
         report(atrow=baserow  ,atcol=j+1,align=center) %strval(100*frac(2),"##.#")
         report(atrow=baserow+1,atcol=j+1,align=center) %strval(100*frac(1),"##.#")+","+%strval(100*frac(3),"##.#")
      end do j
      compute baserow=baserow+2
   end do k
   report(action=show)
end do i
smpl
end

