*
* This is a collection of procedures for setting up a Gibbs sampler for
* a (linear) SUR, such as a near-VAR. @SURGibbsSetup is done outside the
* simulation loop, while @SURGibbsDataInfo and @SURGibbsSigma are used
* inside it. An example and a further explanation are provided in the
* "Bayesian Econometrics" workbook.
*
* @SURGibbsSetup model start end
* examines the model, figures out the positions in the grand B matrix
* for the coefficients and computes the required cross product matrix of
* data over the range <<start>> to <<end>>.
*
* @SURGibbsSigma(bdraw)
* is a function which uses the input draw for the coefficients (<<bdraw>>) and
* computes T x the covariance matrix of residuals at that set of coefficients.
*
* @SURGibbsDataInfo hdraw(input) hdata(output) hbdata(output)
* uses the current value of <<hdraw>> (draw for the residual precision
* matrix) to compute the data evidence on the precision of beta
* (returned as <<hdata>>) and that precision times the mean (returned as
* <<hbdata>>)
*
* The calculations inside the loop will go approximately like this, where
*
* <<hxxprior>> is the prior residual precision x its degrees of freedom
* <<wishdof>> is the posterior degrees of freedom for the residual precision
* (number of observations + degrees of freedom of the prior).
* <<hprior>> is the prior precision for the coefficients
* <<bprior>> is the prior mean for the coefficients
*
*   compute covmat=SURGibbsSigma(bdraw)+hxxprior
*   compute hdraw=%ranwishartf(%decomp(inv(covmat)),wishdof)
*   *
*   @SURGibbsDataInfo hdraw hdata hbdata
*   *
*   compute bdraw=%ranmvposthb(hdata+hprior,hbdata+hprior*bprior)
*
* where HDRAW will be the draw for the precision (inverse
* covariance matrix) of the residuals and BDRAW will be the stacked
* coefficient vector for the model.
*
* Revision Schedule:
*  05/2009 Written by Tom Doan, Estima
*  01/2010 Corrected problem with lookup vector in SURGibbsSetup.
*  09/2018 Change comments to reflect more efficient %RANMVPOSTHB function*

procedure SURGibbsSetup model start end
type model model
type integer start end
*
declare integer           neqn          ;* Number of equations
declare integer           ntotal        ;* Number of coefficients in model
declare vect[int]         esize         ;* Size of each equation
declare vect[int]         cxxdepvar     ;* Positions of dep var's in cross product
declare vect[int]         eqn           ;* Equation of each coefficient
declare vect[int]         eqnpos        ;* Position within equation of coefficient
declare vect[int]         cxxlookup     ;* Position in cross product of each coefficient
declare rect              bblock        ;* [-B/I] array used for computing covariance matrix
*
local rect[int] table(2,0)
local vect[int] dvlist
local integer   i j ipos
local rect[int] etable
local vect[vector[int]] lookup
*
compute neqn=%modelsize(model)
*
dim lookup(neqn)
dim esize(neqn)
*
* Create the full list of explanatory variables without repetitions
* Also, create the lookup list for positions in the list for
* the coefficients in each equation.
*
dim table(2,0)
compute dvlist=%rlempty()
compute ntotal=0
do i=1,neqn
   compute esize(i)=%eqnsize(%modeleqn(model,i))
   compute ntotal=ntotal+esize(i)
   compute etable=%eqntable(%modeleqn(model,i))
   compute table=%tablemerge(table,etable)
   dim lookup(i)(esize(i))
   ewise lookup(i)(j)=%tablefind(table,etable(1,j),etable(2,j))
   compute dvlist=%rladdone(dvlist,%modeldepvars(model)(i))
end do i
*
* Add the dependent variables
*
compute table=%tablemerge(table,%tablefromrl(dvlist))
*
* Compute the grand cross product matrix
*
cmom start end
# %rlfromtable(table)
*
* Set up the BBLOCK matrix. This starts with all zeros, then
* fills in 1.0's for the positions of the dependent variables.
* Also, set up the eqn and cxxlookup arrays for the coefficients.
* eqn(j) is the equation which has coefficient j
* cxxlookup(j) is the position of coefficient j in the cross product
*
compute bblock=%zeros(%ncmom,neqn)
dim eqn(ntotal) eqnpos(ntotal) cxxlookup(ntotal) cxxdepvar(neqn)
compute ipos=0
do i=1,neqn
   compute cxxdepvar(i)=%tablefind(table,dvlist(i),0)
   compute bblock(cxxdepvar(i),i)=1.0
   do j=1,esize(i)
      compute eqn(ipos+j)=i,eqnpos(ipos+j)=j,cxxlookup(ipos+j)=lookup(i)(j)
   end do j
   compute ipos=ipos+esize(i)
end do i
end
***************************************************************
function SURGibbsSigma bdraw
type symm SURGibbsSigma
type vect bdraw
*
local integer ipos i j
*
* Poke (-the) coefficients from bdraw into the proper slots in
* bblock
*
compute ipos=0
do i=1,neqn
   do j=1,esize(i)
      compute bblock(cxxlookup(ipos+j),i)=-bdraw(ipos+j)
   end do j
   compute ipos=ipos+esize(i)
end do i
*
* Compute the quadratic form -> (unscaled) sigma
*
compute SURGibbsSigma=%mqform(%cmom,bblock)
end
*****************************************************************
procedure SURGibbsDataInfo hdraw hdata hbdata
type symm hdraw
type symm *hdata
type vect *hbdata
*
local integer i j
*
dim hdata(ntotal,ntotal) hbdata(ntotal)
*
ewise hdata(i,j)=hdraw(eqn(i),eqn(j))*%cmom(cxxlookup(i),cxxlookup(j))
do i=1,ntotal
   sstats 1 neqn hdraw(eqn(i),t)*%cmom(cxxlookup(i),cxxdepvar(t))>>hbdata(i)
end do i
end
