*
* @ShortAndLong(options) sigma
*
* ShortAndLong computes a factorization of sigma with a combination of
* short and long run restrictions. It can only be applied to
* just-identified parameter- izations where the restrictions are zero
* restrictions. If you have an underidentified parameterization, you can
* compute the "RPERP" matrix which maps free parameters into the loading
* matrix, but you have to use NOESTIMATE.
*
* Parameters:
*  sigma     (input) the covariance matrix of the VAR residuals
*
* The basic model is u = Fv, where v are the structural shocks,
* restricted to be orthonormal, that is, Evv'=I. FF' is thus equal to
* sigma for a just identified model.
*
* Options (inputs):
*  SR=NxN matrix with 0's in the slots in F restricted to be zero.
*  LR=NxN matrix with 0's in the slots in C(1)F restricted to be zero.
*  MASUMS=NxN matrix of the sum of the MAR coefficients of the VAR (C(1))
*
* Options (outputs):
*  FACTOR=factor matrix (F)
*  RPERP=mapping from free parameters to vec(F)
*
*  [ESTIMATE]/NOESTIMATE
*  ITERATIONS=maximum number of iterations [100]
*  CVCRIT=convergence criterion [.000001]. This is checked against the
*    difference between the outer product of the factor and the input
*    covariance matrix.
*  INITIAL=initial guess for the factor matrix F. By default, this takes
*    a Cholesky factor of sigma and picks the closest fit to that.
*  TRACE/NOTRACE
*
*  Example (for a 5 variable system)
*  dec rect lr(5,5) sr(5,5)
*  input lr
*   . . . . .
*   . . 0 0 0
*   . . . 0 0
*   . . . . .
*   . . . . .
*   input sr
*   . 0 0 0 0
*   . . . . .
*   . . . . .
*   . . . . 0
*   . . . . .
*   @ShortAndLong(lr=lr,sr=sr,masum=inv(%varlagsums),factor=f) %sigma
*
* Revision Schedule:
*   12/2005 Written by Tom Doan, Estima
*   07/2007 ESTIMATE, RPERP and FACTOR options added
*   08/2012 Switch to use of %GINV to handle situation with all long-run
*           restrictions.
*   12/2015 Add TRACE option. Restrict size of step for first few iterations.
*   08/2016 Add check for problem with redundant restrictions due to reduced
*           rank MASUMS
*******************************************************************************
procedure ShortAndLong sigma f
type symm sigma
type rect *f
*
option rect    sr
option rect    lr
option rect    masums

option rect    *factor
option rect    *rperp
*
option switch  estimate     1
option integer iterations  100
option real    cvcrit      .000001
option rect    initial
option switch  trace       0
*
local integer  n nr nfree fill i j iters
local rect     r b rp df temp
local vect     theta delta
local symm     gap b2b
local vect[int] rrwzcount
*
compute n=%rows(sigma)

if .not.%defined(sigma) {
   disp "Syntax Error: @ShortAndLong(options) >>>SIGMA<<<"
   return
}
if %defined(sr).and.(%rows(sr)<>n.or.%cols(sr)<>n) {
   disp "@ShortAndLong: SR matrix must be " n "x" n
   return
}
if %defined(lr).and.(%rows(lr)<>n.or.%cols(lr)<>n) {
   disp "@ShortAndLong: LR matrix must be " n "x" n
   return
}
if %defined(lr).and..not.%defined(masums) {
   disp "@ShortAndLong: LR matrix needs MASUMS option"
   return
}
if %defined(masums).and.(%rows(masums)<>n.or.%cols(masums)<>n) {
   disp "@ShortAndLong: MASUMS matrix must be " n "x" n
   return
   }
*
* Compute the number of restrictions. Also, count the number of
* restrictions in each column to do the RRWZ test
*
compute nr=0
dim rrwzcount(n)
ewise rrwzcount(i)=0
do i=1,n
   do j=1,n
      if %defined(sr).and.sr(i,j)==0
         compute nr=nr+1,rrwzcount(j)=rrwzcount(j)+1
      if %defined(lr).and.lr(i,j)==0
         compute nr=nr+1,rrwzcount(j)=rrwzcount(j)+1
   end do j
end do i
*
compute nfree=n^2-nr
compute r=%zeros(nr,n^2)
compute fill=0
do i=1,n
   do j=1,n
      if %defined(sr).and.sr(i,j)==0 {
         compute fill=fill+1
         compute r(fill,(j-1)*n+i)=1.0
      }
      if %defined(lr).and.lr(i,j)==0 {
         compute fill=fill+1
         compute %psubmat(r,fill,(j-1)*n+1,%xrow(masums,i))
      }
   end do j
end do i
compute rp=%perp(r)
if %defined(rperp)
   compute rperp=rp
if .not.estimate
   return
*
* If we want to estimate using this procedure, we need a just-identified
* model
*
if nr<>n*(n-1)/2 {
   disp "##@ShortAndLong Can Only Estimate Just-Identified Models"
   disp "  Had" fill "Restrictions. Need" n*(n-1)/2
   return
}
*
* Check whether we had redundant restrictions (R was deficient rank)
*
if %cols(rp)>nfree {
   disp "##@ShortAndLong inputs include " %cols(rp)-nfree "redundant restriction(s)"
   disp "  Typically this is due to LR restrictions which aren't separate restrictions"
   disp "  because of a reduced rank MASUMS matrix (in a VECM)"
   return
}
*
* Check and warn on the RRWZ conditions.
*
do i=0,n-1
   do j=1,n
      if rrwzcount(j)==i
         goto rrwzok
   end do j
   disp "##@ShortAndLong Rubio-Ramirez-Waggoner-Zha condition fails. Model isn't globally identified."
   disp "  Have Columns with " rrwzcount "zeros. Need " %seq(0,n-1) "in some order"
   return
:rrwzok
end do i
*
* theta is the vector of free coefficients. It's initialized to the
* closest fit to the "b" vector which is the vectorized Cholesky factor
* of sigma (unless the user inputs a different initial matrix.)
*
dim theta(nfree)

if %defined(initial)
   compute b=initial
else
   compute b=%vec(%decomp(sigma))
compute theta=inv(tr(rp)*rp)*tr(rp)*b
*
dim df(nfree,nfree)
*
* Solves the system of equations sigma-b(theta)*tr(b(theta))=0 by
* Newton's method
*
do iters=1,iterations
   compute b=%vectorect(rp*theta,n)
   *
   compute gap=b*tr(b)-sigma
   compute %cvcrit=sqrt(%normsqr(gap))
   if trace
      disp "Iteration" iters "cvcrit" %cvcrit
   if %cvcrit<cvcrit
      break
   *
   * Computes the matrix of derivatives of the "nfree" components of the
   * VEC'ed objective function by the components of theta.
   *
   do i=1,nfree
      compute temp=b*tr(%vectorect(%xcol(rp,i),n))
      compute [symm] b2b=temp+tr(temp)
      compute %psubmat(df,1,i,%vec(b2b))
   end do i
   *
   * Compute the Newton's method step, and adjust theta
   *
   compute delta=%ginv(df)*%vec(gap)
   *
   * Take a partial step on the first few iterations
   *
   if iters<=4
      compute theta=theta-(1-exp(-.10*iters*iters))*delta
   else
      compute theta=theta-delta
end do iters
compute %iters=iters
compute %converged=fix((%cvcrit<cvcrit))
if .not.%converged
   disp "ShortAndLong failed to converge to a solution"
if %defined(f)
   compute f=b
if %defined(factor)
   compute factor=b
end shortandlongrun
