*
* CONDITION.SRC
* Computes conditional forecasts. See Section 7.7.2 of the RATS User's Guide.
* https://estima.com/docs/RATS%2010%20Users%20Guide.pdf#page=260
*
*  Syntax:
*
*  @CONDITION( options ) # of constraints
*  #  series   entry    value       (one for each constraint)
*
*     or
*
*  @CONDITION(GENERAL,other options )  # of constraints
*  #  vector with linear combination of endogenous variables (one pair
*  #  entry   value                                           for each constraint)
*
*  Options:
*    model       = VAR model used for forecasts
*    from        = start period for forecasts [last estimation + 1]
*    to          = end period for forecasts [until last constrained value]
*    steps       = number of forecast steps (alternative to using TO option)
*    cv          = covariance matrix of residuals [model's]
*    factor      = factor of covariance matrix [not used]
*    simulate/[nosimulate] if simulated, draws randomly from the constrained density
*
*    results     = VECT[SERIES] containing the constrained forecasts (output)
*    shocks      = VECT[SERIES] with the orthogonalized shocks used (output)
*
*  Supplementary Cards (first form, one for each constraint)
*    series      = endogenous variable constrained
*    entry       = entry of constraint
*    value       = target value
*
*  Supplementary Cards (second form (GENERAL option), 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
*
*  Revision Schedule
*    12/2003 SIMULATE option added
*    01/2005 GENERAL option added
*    07/2007 Rewritten to use v7 syntax
*    11/2012 Add test for constraint < forecast start
*    04/2019 Adapt to allow FACTOR matrix with fewer columns
*
procedure condition nconstr
type integer 	nconstr
*
option switch  general      0
option switch  simulate     0
option model   model
option integer	steps
option integer from
option integer to
option symm    cv
option rect    factor
option vect[series] *results
option vect[series] *shocks
*
local integer  nvar nshocks nbeg nfore
local index    step var
local vector   value
local integer 	i j nsteps nfree constr consentr
local integer 	fails
local rect     capr uover vover vfactor overr
local vect     lowr u v
local rect[series] impfunc
local vect[vect] lcomb
dec symm 		   rrrr
*
* Process supplementary cards for PROCEDURE.
*
dimension step(nconstr) var(nconstr) value(nconstr) lcomb(nconstr)
*
if .not.%defined(model) {
   disp "Error: @Condition(MODEL=MODEL TO FORECAST,other options) # of constraints"
   return
}
*
if .not.%defined(from)
   compute nbeg=%regend()+1
else
   compute nbeg=from
*
compute nvar=%modelsize(model)
compute nsteps=0
compute fails =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)<1 {
      disp "Error: @CONDITION. Constraint" i "is before forecast start period" %datelabel(nbeg)
      compute fails=1
   }
   if step(i)>nsteps  ; compute nsteps=step(i)
end do i
if fails
   return
*
if .not.%defined(steps).and..not.%defined(to)
   compute nfore=nsteps
else
if %defined(to)
   compute nfore=to-nbeg+1
else
   compute nfore=steps
*
* Replace VAR(I) with the equation which has VAR(I) as its dependent variable.
*
if .not.general {
   do i=1,nconstr
      compute var(i)=%modelfind(model,var(i))
      if var(i)==0 {
         display "Constraints Must Use One of the Endogenous Variables"
         display "Please Check Constraint Number" i
         return
      }
   end do i
}
*
if %defined(factor)
   compute vfactor=factor
else
if %defined(cv)
   compute vfactor=%decomp(cv)
else
   compute vfactor=%decomp(%modelgetvcv(model))
*
compute nshocks=%cols(vfactor)
forecast(model=model,results=results,steps=nsteps,from=nbeg)
impulse(model=model,results=impfunc,noprint,factor=vfactor,steps=nsteps)
*
compute nfree=nshocks*nsteps
dimension capr(nfree,nconstr)
dimension u(nvar*nsteps) v(nfree) lowr(nconstr)
do constr=1,nconstr
   if .not.general
      compute lowr(constr)=value(constr)-results(var(constr))(step(constr)+nbeg-1)
   else
      compute lowr(constr)=value(constr)-%dot(lcomb(constr),%xt(results,step(constr)+nbeg-1))
   overlay capr(1,constr) with overr(nsteps,nshocks)
   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,nshocks)
*
* If simulations are requested, do a draw from the constrained density for the
* orthogonal shocks.
*
if simulate {
   compute rrrr=%identity(nsteps*nshocks)-capr*inv(tr(capr)*capr)*tr(capr)
   compute v=v+%ranmvnormal(%decomp(rrrr))
}
*
* Transform to non-orthogonal U
*
compute uover=vover*tr(vfactor)
*
* Save the orthogonalized shocks, if requested
*
if %defined(shocks) {
   dim shocks(nshocks)
   do i=1,nshocks
      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=results,steps=nfore,from=nbeg)
end condition
