*
*  CONDITION.SRC
*  Computes conditional forecasts. See Section 10.14 of the Manual
*
*  Syntax:
*
*  @CONDITION( options )  model  steps  start  VCV Matrix  constraints  forecasts
*  #  series   entry    value  (one for each constraint)
*
*     or
*
*  @CONDITION(GENERAL,other options )  model  steps  start  VCV Matrix  constraints  forecasts
*  #  vector with linear combination of endogenous variables (one pair for each constraint)
*  #  entry value
*
*    model       = VAR model used for forecasts
*    steps       = number of forecast steps
*    start       = start period of forecast
*    VCV Matrix  = SYMMETRIC covariance matrix of residuals. The rows must
*                  correspond to equations in order
*    constraints = number of constraints on future values
*    forecasts   = VECT[SERIES] containing the constrained forecasts. This
*                  is returned by the procedure.
*
*  Supplementary Cards, first type (one for each constraint)
*    series      = endogenous variable constrained
*    entry       = entry of constraint
*    value       = target value
*
*  Supplementary Cards, second type (one of each for each constraint)
*    vector ...  = VECTOR or real numbers/expressions with linear combination of
*                    endogenous variables to constrain
*    entry       = entry of constraint
*    value       = target value for the linear combination
*
*  Options
*    general/[nogeneral] use GENERAL to use constraints on linear combinations of series, rather
*                        than on individual variables
*    shocks = VECT[SERIES] with the orthogonalized shocks used
*    decomp = factorization matrix [Cholesky factor of VCV matrix]
*    simulate/[nosimulate] if simulated, draws randomly from the constrained density
*
*  Revision Schedule
*    12/03 SIMULATE option added
*    01/05 GENERAL option added
*    03/08 Renamed CONDITION from CNDITION
*
procedure condition  model nfore nbeg vmat nconstr fcstout
type model model
type integer nfore nbeg nconstr
type vect[series] *fcstout
type symmetric vmat
*
option switch        simulate 0
option rect          decomp
option vect[series] *shocks
option switch        general 0
*
local integer nvar
local index step var
local vector value
local integer i j nsteps ieqn iblock nfree constr consentr enumber
local rect capr uover vover vdecomp overr fcsts impulses
local vect lowr u v
local rect[series] impfunc
local vect[integer] depvars
local vect[vect] lcomb
local symm rrrr
*
*   Process  supplementary cards for PROCEDURE.
*
dimension step(nconstr) var(nconstr) value(nconstr) lcomb(nconstr)
*
compute nvar=%modelsize(model)
compute nsteps=0
do i=1,nconstr
  if .not.general
     enter var(i) step(i) value(i)
  else {
     dim lcomb(i)(nvar)
     enter lcomb(i)
     enter step(i) value(i)
  }
  compute step(i)=step(i)-nbeg+1
  if step(i)>nsteps  ; compute nsteps=step(i)
end do i
*
*  Replace VAR(I) with the equation which has VAR(I) as its dependent
*  variable.
*
compute depvars=%modeldepvars(model)
if .not.general {
   do i=1,nconstr
      compute enumber=0
      do j=1,nvar
         if var(i)==depvars(j) {
            compute enumber=j
            break
         }
      end do j
      if enumber==0 {
         display 'Constraints Must Use One of the Endogenous Variables'
         display 'Please Check Constraint Number' I
         return
      }
      compute var(i)=enumber
   end do i
}

if %defined(decomp)
   compute vdecomp=decomp
else
   compute vdecomp=%decomp(vmat)

forecast(model=model,results=fcstout) * nsteps nbeg
impulse(model=model,results=impfunc,noprint,decomp=vdecomp) * nsteps *

compute nfree=nvar*nsteps
dimension capr(nfree,nconstr) vdecomp(nvar,nvar)
dimension u(nfree) v(nfree) lowr(nconstr)
do constr=1,nconstr
   if .not.general
      compute lowr(constr)=value(constr)-fcstout(var(constr))(step(constr)+nbeg-1)
   else
      compute lowr(constr)=value(constr)-%dot(lcomb(constr),%xt(fcstout,step(constr)+nbeg-1))
   overlay capr(1,constr) with overr(nsteps,nvar)
   compute consentr=1+step(constr)
   if .not.general
      ewise overr(i,j)=%if(i<consentr,impfunc(var(constr),j)(consentr-i),0.0)
   else
      ewise overr(i,j)=%if(i<consentr,%dot(lcomb(constr),%xcol(%xt(impfunc,consentr-i),j)),0.0)
end do constr
*
*  Compute optimal orthogonalized V
*
compute v=capr*inv( tr(capr)*capr )*lowr
overlay u(1) with uover(nsteps,nvar)
overlay v(1) with vover(nsteps,nvar)
*
*  If simulations are requested, do a draw from the constrained
*  density for the orthogonal shocks
*
if simulate {
   compute rrrr=%identity(nsteps*nvar)-capr*inv(tr(capr)*capr)*tr(capr)
   compute v=v+%ranmvnormal(%decomp(rrrr))
}
*
*  Transform to non-orthogonal U
*
compute uover=vover*tr(vdecomp)
*
*  Save the orthogonalized shocks
*
if %defined(shocks) {
   dim shocks(nvar)
   do i=1,nvar
      set shocks(i) nbeg nbeg+nsteps-1 = vover(t-nbeg+1,i)
   end do i
}
*
*   Forecast with U shocks
*
forecast(model=model,matrix=uover,results=fcstout) * nfore nbeg
end condition

