*
* @MCProcessIRF(model=model used in generating responses,other options )
*
* Extracts information for graphing or printing tables for 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.
*
* NVAR is either taken from the model input to the procedure or by using
* the NVAR option.
*
* Options:
*
* MODEL=model used in generating responses
* or
* NVAR=number of response variables
*
* CENTER=[MEAN]/MEDIAN/INPUT
* 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.
*
* WEIGHTS=series of (relative) weights on the draws (from importance
* sampling) [equal weights]. Note that this can dramatically slow down
* the calculation if you ask for percentiles and have many (5000+)
* draws.
*
* ACCUMULATE=list of variables (by original position) to accumulate [none]
* If you already did accumulation when you generated the impulse responses,
* don't repeat it here.
*
* These are the outputs. They all need to be used. Each generates a
* RECT[SERIES], where x(i,j) is the series of responses of variable i to
* shock j.
*
* IRF=RECT[SERIES] with the central measure of the impulse responses
* LOWER=RECT[SERIES] with the lower bounds
* UPPER=RECT[SERIES] with the upper bounds
*
* Revision Schedule:
* 11/2009 Stripped from MCGraphIRF.src procedure
* 03/2011 Add WEIGHTS option
* 07/2013 Change to allow NVAR option to be used rather than MODEL.
* 10/2016 Add ACCUM option.
*
procedure MCProcessIRF
*
option model model
option vector percentiles ||.16,.84||
option real stddev
option choice center 1 mean median input
option rect[series] impulses
option rect[series] *upper
option rect[series] *lower
option rect[series] *irf
option series weights
option vect[integer] accum
option integer nvar
*
local integer draws steps nshocks npercent
local real sigma worksum totalweights
local integer i j k ia kk
local vector work
local vector frac
local integer ntargets
local vect request
*
declare vect[rect] %%responses
*
if %rows(%%responses)==0 {
disp "####@MCProcessIRF - needs procedure to draw responses first (%%responses not defined)"
return
}
if .not.%defined(model).and..not.%defined(nvar) {
disp "Syntax: @MCProcessIRF needs either MODEL option or NVAR option"
return
}
if center==3.and..not.%defined(impulses) {
disp "####@MCProcessIRF(CENTER=INPUT) needs IMPULSES option"
return
}
if .not.%defined(lower) {
disp "####@MCProcessIRF(LOWER=RECT[SERIES] for lower bounds,other options)"
return
}
if .not.%defined(upper) {
disp "####@MCProcessIRF(UPPER=RECT[SERIES] for upper bounds,other options)"
return
}
if .not.%defined(irf) {
disp "####@MCProcessIRF(IRF=RECT[SERIES] for central values,other options)"
return
}
*
* Take the number of variables (targets) out of the model
*
if %defined(model)
compute ntargets=%modelsize(model)
else
compute ntargets=nvar
*
* 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)
*
* 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=%size(percentiles)
if center==2 {
if npercent>0
compute request=percentiles~~||.50||
else
compute request=||.50||
}
else
if npercent>0
compute request=percentiles
dim upper(ntargets,nshocks) lower(ntargets,nshocks) irf(ntargets,nshocks)
do i=1,ntargets
do j=1,nshocks
set irf(i,j) 1 steps = 0.0
set lower(i,j) 1 steps = 0.0
set upper(i,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+i,k)
*
* See if we're supposed to accumulate the response. If so, add in
* the responses of the preceding steps (for each replication).
*
if %defined(accum) {
do ia=1,%size(accum)
if i==accum(ia) {
do kk=1,k-1
ewise work(t)=work(t)+%%responses(t)((j-1)*ntargets+i,kk)
end do kk
break
}
end do ia
}
*
* Compute whatever percentiles are needed
*
if %defined(weights)
compute frac=%wfractiles(work,weights,request)
else
compute frac=%fractiles(work,request)
*
* Choose the central value
*
if center==1.and.%defined(weights) {
sstats 1 draws work(t)*weights>>worksum weights>>totalweights
compute irf(i,j)(k)=worksum/totalweights
}
else
if center==1
sstats(mean) 1 draws work(t)>>irf(i,j)(k)
else
if center==2
compute irf(i,j)(k)=frac(npercent+1)
else
compute irf(i,j)(k)=impulses(i,j)(k)
if %defined(stddev) {
if %defined(weights) {
sstats 1 draws (work(t)-irf(i,j)(k))^2*weights>>worksum
compute sigma=sqrt(worksum/totalweights)
}
else {
sstats(mean) 1 draws (work(t)-irf(i,j)(k))^2>>worksum
compute sigma=sqrt(worksum)
}
compute lower(i,j)(k)=irf(i,j)(k)-stddev*sigma
compute upper(i,j)(k)=irf(i,j)(k)+stddev*sigma
}
else {
compute lower(i,j)(k)=frac(1)
compute upper(i,j)(k)=frac(2)
}
end do k
end do j
end do i
*
end MCProcessIRF