*
* @VARIRFDelta( options ) covxx
* computes the covariance matrix of the IRF to a fixed shock using the
* delta method.
*
* Parameters:
*   covxx   (required) covariance matrix of the responses at the indicated horizon.
*
* Options:
*   model=VAR model (required)
*   h    =horizon desired (one per call) [1]
*   shock=VECTOR with initial shock [unit shock to first component]
*
* Revision Schedule:
*  11/2009 Written by Tom Doan, Estima
*  11/2018 Add i j as local variables
*  07/2022 Correct problem with handling of unit shock
*
procedure VARIRFDelta covxx
type symm *covxx
*
option model   model
option vector  shock
option integer h       1
*
local rect       a
local symm       covxx
local rect       derive basedd
local rect[vect] dgrand
local integer    nreg nvar nlags
local integer    i j k l hh
local integer    i1 i2 j1 j2
local vect[int]  reorder
local vector     shockw
*
if .not.%defined(model).or..not.%defined(covxx) {
   disp "@VARIRFDelta(model=VAR model,other options)   covxx"
   return
}
compute a=%modelcompanion(model)
*
compute nvar =%modelsize(model)
compute nreg =%rows(a)
compute nlags=nreg/nvar
*
* Shuffles the components of the companion matrix to the order of
* estimation.
*
dim reorder(nreg)
ewise reorder(i)=(%clock(i,nvar)-1)*nlags+(i-1)/nvar+1
*
dim covxx(nvar,nvar)
dim derive(nreg,nreg) basedd(nreg,nreg)
dim dgrand(nvar,nreg)
*
* Use a unit vector if the user hasn't provided an alternative.
*
if %defined(shock)
   compute shockw=shock
else
   compute shockw=%unitv(nvar,1)

do i=1,nvar
   do j=1,nreg
      ewise basedd(k,l)=(k==i.and.l==j)
      compute derive=%zeros(nreg,nreg)
      do hh=1,h
         compute derive=derive*a+a^(hh-1)*basedd
      end do hh
      compute dgrand(i,j)=%xsubmat(derive,1,nvar,1,nvar)*shockw
   end do j
end do i
*
compute covxx=%zeros(nvar,nvar)
do i1=1,nvar
   do i2=1,nvar
      do j1=1,nreg
         do j2=1,nreg
            compute covxx=covxx+%sigma(i1,i2)*%xx(reorder(j1),reorder(j2))*$
               dgrand(i1,j1)*tr(dgrand(i2,j2))
         end do j2
      end do j1
   end do i2
end do i1
end
