*
* @BVarBuildPriorMN(options)  hbpriors hpriors
* builds a "Minnesota" prior for a VAR model.
*
* Parameters (required)
*   HBPRIORS   VECT[VECT] accumulation of precision matrix x prior mean
*              vectors for the elements of the prior for each equation.
*   HPRIORS    VECT[SYMM] accumulation of the precision matrices for the
*              elements of the prior for each equation.
*
* Options:
*   MODEL=MODEL for which the prior is being constructed
*
*   SCALEFACTORS=VECTOR of OLS autoregression scale factors
*   [CALCULATE]/NOCALC
*      If you need to set up several priors based upon a single set of data,
*      you can use NOCALC to avoid doing the calculation of the cross product
*      matrix and the scale factors. You need to have the scale factors the
*      first time through. Then use the same option on the future calls.
*
*   EQUATION=equation defining the form of the estimated equation
*   INIT/[NOINIT] initialize (dimension and zero) the hbprior and hprior
*
*   TIGHT=overall tightness parameter [.20]
*   OTHER=relative tightness on other variables [1.0]
*   MATRIX=RECTANGULAR array of weights for general handling of variables.
*   CONTIGHT=relative tightness on deterministic variables [flat]
*   DECAY=harmonic decay exponent (std dev is l^-decay) [0=no decay]
*   MEAN=prior mean for first lags [1.0]
*   MVECTOR=VECTOR of prior means (overrides MEAN if used) of individual equations.
*
*  Revision Schedule:
*   05/2009 Written by Tom Doan, Estima
*   08/2013 Add MATRIX option
*   04/2017 Correct the lag decay (was reciprocal of what it should be)
*   03/2024 General clean up to make structure of input model more general
*           (allows skipped lags, not fully identical equations).
*
procedure BVARBuildPriorMN hbpriors hpriors
type vect[vect] *hbpriors
type vect[symm] *hpriors
*
option model  model
*
option vector *scalefactors
option switch calculate 1
*
option real   tight    .20
option real   other    1.0
option real   contight
option rect   matrix
option real   decay    0.0
option real   mean     1.0
option vector mvector
*
local  table   table
local  integer nvar ncoeff nlags sdet lag i j k
local  real    meanb se
local  vect    lagmeans olssee
local  real    fij
local  equation eqn
local  vect[int] laglist detlist reglist
*
if .not.%defined(model) {
   disp "###@BVARBuildPriorMN(MODEL=var model,...)"
   return
}
compute eqn=%modeleqn(model,1)
compute table=%eqntable(eqn)
compute nvar=%modelsize(model)
if %defined(matrix) {
   if %rows(matrix)<>nvar.or.%cols(matrix)<>nvar {
      disp "###@BVARBuildPriorMN(MATRIX option) error"
      disp "MATRIX is " %rows(matrix) "x" %cols(matrix) ". Must be" nvar "x" nvar
      return
   }
}
*
* Figure out the highest lag of the first dependent variable. At the
* same time, extract the "deterministic" variables.
*
dim detlist(0)
compute ncoeff=%cols(table)
compute nlags=0
do i=1,ncoeff
   if table(1,i)==%eqndepvar(eqn)
      compute nlags=%imax(nlags,table(2,i))
   if %modelfind(model,table(1,i)==0)
      compute detlist=%rladdlag(detlist,table(1,i),table(2,i))
end do i
*
compute sdet=nvar*nlags+1
dim olssee(nvar)
*
* If needed, compute the regression standard deviations for univariate autoregressions
*
if calculate {
   do i=1,nvar
      compute eqn=%modeleqn(model,i)
      compute laglist=%tablefindlags(%eqntable(eqn),%eqndepvar(eqn))
      compute reglist=%rladdlaglist(detlist,%eqndepvar(eqn),laglist)
      linreg(noprint) %eqndepvar(eqn)
      # reglist
      compute olssee(i)=sqrt(%seesq)
   end do i
   if %defined(scalefactors)
      compute scalefactors=olssee
}
else {
   if .not.%defined(scalefactors) {
      disp "###@BVARBuildPriorMN(NOCALC,SCALEFACTOR=VECTOR of scale factors,...)"
      return
   }
   if %size(scalefactors)<>nvar {
      disp "###@BVARBuildPriorMN(SCALEFACTOR=...) has wrong dimensions"
      return
   }
   compute olssee=scalefactors
}
*
dim hbpriors(nvar) hpriors(nvar)
*
if %defined(mvector)
   compute lagmeans=mvector
else
   compute lagmeans=%fill(nvar,1,mean)
*
* Prior for equation i
*
do i=1,nvar
   *
   * Initialize the prior matrices
   *
   compute eqn=%modeleqn(model,i)
   compute table=%eqntable(eqn)
   @BVARBuildPrior(init,equation=eqn) hbpriors(i) hpriors(i)
   *
   do k=1,%eqnsize(eqn)
      compute j=%modelfind(model,table(1,k))
      *
      * If this is *not* a dependent variable of any equation, apply the
      * "contight" tightness (0=infinite, which means that we skip it since
      * zero precision is the default).
      *
      if j==0 {
         if %defined(contight)
            @BVARBuildPrior(position=k,mean=0.0,stderr=tight*contight) hbpriors(i) hpriors(i)
      }
      else {
         *
         *  See if this is the 1st lag of the dependent variable of equation i.
         *  If so, the mean is the request 1st lag mean. If not, the mean is zero.
         *
         *
         compute lag=table(2,k)
         if j==i.and.lag==1
            compute meanb=lagmeans(i)
         else
            compute meanb=0.0
         if %defined(matrix)
            compute fij=matrix(i,j)
         else
            compute fij=%if(i==j,tight,tight*other)
         compute se=fij*(olssee(i)/olssee(j))*(lag^(-decay))
         @BVARBuildPrior(position=k,mean=meanb,stderr=se) hbpriors(i) hpriors(i)
      }
   end do k
end do i
end
*********************************************************
*
* Converts the output from @BVARBuildPriorMN into a stacked mean and blocked
* precision matrix.
*
procedure BVARFinishPrior hbpriors hpriors bprior hprior
type vect[vect] hbpriors
type vect[symm] hpriors
type rect       *bprior
type symm       *hprior
*
local integer i nvar
compute nvar=%rows(hbpriors)
local vect[rect] hpriorsr(nvar)
local rect[vect] bpriorsr(1,nvar)
do i=1,nvar
   compute hpriorsr(i)=hpriors(i)
   compute bpriorsr(1,i)=%ginv(hpriors(i))*hbpriors(i)
end do i
compute hprior=%blockdiag(hpriorsr)
compute bprior=%blockglue(bpriorsr)
end
