The RATS Software Forum

Code: Select all

compute meanx2=%digamma(0.5)-log(0.5)
compute varx2 =%trigamma(0.5)
*
*  Create the adjusted log squared return series
*
dec vect[series] w(2)
set dlogp = log(price/price{1})
diff(center) dlogp / demean
set w(1) = log(demean**2)-meanx2
set dlogp = log(ADMD/ADMD{1})
diff(center) dlogp / demean
set w(2) = log(demean**2)-meanx2

*
dec packed rhosd(1,1) sigetaf(2,2)
dec symm   sv(2,2)
*
* %rhosw combines the packed lower triangular "rho" matrix with the
* known variances of the components to give the full covariance matrix.
* The sigma eta matrix is easier to estimate in square root form so it
* is also a packed lower triangle. This also makes the estimation
* process more similar to the factor method later.
*
function %rhosw
type symm %rhosw
local integer i j
dim %rhosw(2,2)
ewise %rhosw(i,j)=%if(i==j,varx2,varx2*rhosd(i-1,j))
end
*
function %ltouterxxF
type symm %ltouterxxF
local rect lt(2,2)
local integer i j
ewise lt(i,j)=%if(i>=j,sigetaf(i,j),0.0)
compute %ltouterxxF=lt*tr(lt)
end
*
* sigmaeps backs out the covariance matrix of epsilon from the estimate
* for eta
*
function %sigmaeps
type symm %sigmaeps
*
local integer i j n
local real    rhowork rho factor
local parmset rootf
*
dim %sigmaeps(2,2)
do i=1,2
   compute %sigmaeps(i,i)=%pi**2/2.0
   do j=1,i-1
      nonlin(parmset=rootf) rho
      compute rho=rhosd(i-1,j)
      find(noprint,parmset=rootf) root rhosd(i-1,j)-rhowork
         compute rhowork=0.0
         compute factor=2.0*rho**2
         do n=1,100
            compute rhowork=rhowork+factor/n
            compute factor=factor*n/(n+.5)*rho**2
         end do n
         compute rhowork=rhowork*2.0/%pi**2
      end find
      compute %sigmaeps(i,j)=rho*%pi**2/2.0
   end do j
end do i
end
*
* The estimates of the sigmaeta matrix that we get are somewhat
* different from those in the paper, and the likelihood (after adjusting
* for the constants) is higher here. It appears that the maximization
* algorithm used in doing the paper ran up against the invertibility
* boundary (note the .99 in the (2,4) slot in the correlation matrix for
* eta), and either stalled out or hit a slightly lower local maximum.
*
nonlin rhosd sigetaf
compute [rect] c=%identity(2)
compute rhosd=%const(0.0)
compute sigetaf=%sqrt(%diag(||0.001,0.001||))
dlm(y=%xt(w,t),c=c,startup=%(sv=%rhosw(),sw=.001*%ltouterxxF()),sw=sw,sv=sv,$
  presample=diffuse,condition=1,method=bfgs,iters=100) 2 *
*
* Get the covariance matrix for epsilon
*
compute sigmaeps=%sigmaeps()
*
* Show the correlation matrices for epsilon and eta
*
disp %cvtocorr(sigmaeps)
disp %cvtocorr(sw)
*
* Table 4 analysis. Because the estimated sigmaeta is quite a bit
* different, these results are as well. The eigenvectors shown in this
* table are actually scaled by the square roots of eigenvalues as is
* done in table 5, so there is no need for a separate table.
*
eigen(scale) sw eigval eigvect
report(action=define)
report(atrow=1,atcol=1) "(a)"	"Eigenvalues" eigval
report(atrow=2,atcol=2) "Eigenvectors" eigvect
report(atrow=6,atcol=2) "Percentage of Variance" eigval/%sum(eigval)
*
eigen(scale) %cvtocorr(sw) eigval eigvect
report(atrow=8,atcol=1) "(b)"	"Eigenvalues" eigval
report(atrow=9,atcol=2) "Eigenvectors" eigvect
report(atrow=13,atcol=2) "Percentage of Variance" eigval/%sum(eigval)
report(action=show)
*