/*
***
*** Change for v. 5.0 compatibility:
***  In v. 5.0, the VCV instruction sets %SIGMA (even if MATRIX option
***  is used). Unfortunately, CATS uses %SIGMA as a user-defined variable
***  name. With RATS 4., VCV didn't set %SIGMA, so this wasn't a problem,
***  but with RATS 5, the %SIGMA as set by CATS will be replaced by
***  the computed VCV matrix, leading to erroneous results when
***  the RESCHK procedure (Residual Analysis) is called.
***  Short-term fix is to save %SIGMA into %%SIGMA, do the VCV
***  then reset %SIGMA equal to %%SIGMA to restore original values.
***  Long term: recode CATS to use %%SIGMA throughout.
***
***  By Estima, Dec. 1, 2000
***
CATSMISC.SRC  95.02.07
The file is called as a source file by the procedure cats in cats.src
The file contains 5 procedures for CATS.
The procedures are compiled if the option misc is included (default is on)
*/
frequency 0 %DIMZ0*%LAG
*
PROCEDURE %GRAF%
/*
 Produces graphs of the series, residuals and coint. vectors
*/
loc rect ci ciu betark uhat betarku fitmat
loc vect Vhelp uhatmean
loc vect[series] crosscorrel
loc int corrend count i1 i ng selection number
loc string header header1
loc series ser1 ser2 ser3 shad up low actual fitted nval xval yval
loc real win range
loc vect[integer] selectlist
loc vect[string] keylab graphlabel catserlabel alphalabel
*
compute number = 1
dim crosscorrel(%DIMZ0**2) Vhelp(%N) uhatmean(%dimZ0) keylab(2)
smpl %start+%LAG %end
* comp  corrend = %N/4
if %n>20000
 compute corrend = 100
else
 compute corrend = %N/4

set up  1 corrend+1 = 2.0/sqrt(%N1)
set low 1 corrend+1 = -2.0/sqrt(%N1)
comp ci = tr(%beta1)*%ZkMAT
comp ciu= tr(%beta0)*%ZkMAT
comp betaRk = ci - (tr(%beta1)*tr(%M1k)*%iM11)*%Z1MAT
comp betaRku = ciu - (tr(%beta0)*tr(%M1k)*%iM11)*%Z1MAT
comp fitmat = (tr(%M10) - %pimat*tr(%M1k))*%iM11*%Z1MAT + %pimat*%ZkMAT
comp uhat = %Z0MAT - fitmat
do %i=1,%DIMZ0
 ewise Vhelp(i) = uhat(%i,i)
 comp uhatmean(%i) = %sum(Vhelp)*%iN
end do %i
ewise uhat(i,j) = (uhat(i,j)-uhatmean(i))/sqrt(%sigma(i,i))
*
loop
 cntrl(clear)
 *
 dim catserlabel(%rows(%catser)) alphalabel(%dimZ0)
 ewise catserlabel(i) = %catserlab(i)
 ewise alphalabel(i) = %alphalab(i)
 if %runtype == 1 ; dim graphlabel(7)
 else             ; dim graphlabel(6)
 comp graphlabel(1) = 'The Time Series'
 comp graphlabel(2) = 'Residuals'
 comp graphlabel(3) = 'Correlations'
 comp graphlabel(4) = 'The Ci relations'
 comp graphlabel(5) = 'The Eigenvectors'
 comp graphlabel(6) = 'Save the data'
 if %runtype == 1
  if %clock(%gflag,2) == 1
   comp graphlabel(7) = 'Save plots (ON)'
  else ; begin
   comp graphlabel(7) = 'Save plots (OFF)'
   close plot
  end
 *
 select(prompt='Graphics selection',strings=graphlabel,status=%stop) selection
 if %stop==0 ; branch exit
 if selection==1 ; branch ts
 if selection==2 ; branch resids
 if selection==3 ; branch corrs
 if selection==4 ; branch cointvect
 if selection==5 ; branch vvectors
 if selection==6 ; branch savedat
 if selection==7 ; branch saveplot
*
:ts
 select(prompt='Select one or more series',strings=catserlabel,status=%stop) $
 number selectlist
 if %stop == 0 .or. number==0 ; next
 do %i=1,number
  compute %whatg = selectlist(%i)
  set ser1 %start %end = [SERIES]%CATSER(%whatg)
  set ser2 %start+1 %end = ser1(t) - ser1(t-1)
  set ser2 %start %start = %NA
  disp(store=header) %CATSERLAB(%whatg)
  @%GRAFON%
  spgraph(vfield=2)
   graph(header=header,sub='LEVEL') 1
   # ser1
   graph(sub='DIFFERENCE',samesize) 1
   # ser2 / 1
  spgraph(done)
 end do %i
 branch ts
 next
*
:resids
 select(prompt='Select one or more series',strings=alphalabel,status=%stop) $
 number selectlist
 if %stop == 0 .or. number==0 ; next
 do %i=1,number
  compute %whatg = selectlist(%i)
  set actual / = %Z0MAT(%whatg,t-%start-%lag+1)
  set fitted / = FITMAT(%whatg,t-%start-%LAG+1)
  set ser1 / = UHAT(%whatg,t-%start-%LAG+1)
  set ser3 / = ser1
  order ser3 %start+%lag %end
  comp win = 1.06/(%N1**0.2)
  comp range = ser3(%end) - ser3(%start+%lag) + 4*win
  comp ng = fix(range/win)
  set xval 1 1 = ser3(%start+%lag) - win
  set xval 2 ng = xval{1} + win
  set yval 1 ng = 0.0
  comp i1 = %start+%lag
  do i=1,ng
   while ser3(i1) .le. xval(i) .and. i1 .le. %end ; begin
    set yval i i = yval(i) + 1
    comp i1 = i1+1
   end while
  end do i
  set yval 1 ng = yval/(%N1*win)
  set nval 1 ng = %density(xval-0.5*win)
  *
  correlate(noprint,number=corrend) ser1 / ser2
  set shad 1 corrend+1 = ABS(ser2) >=  up(1)
  disp(store=header) 'Actual and Fitted for' %alphalab(%whatg)
  comp keylab(1)  = 'Normal'
  comp keylab(2)  = %alphalab(%whatg)
  @%GRAFON%
  spgraph(vfields=2,hfieds=2)
   graph(header=header) 2
   # actual ; # fitted
   graph(header='Standardized Residuals',samesize) 1
   # ser1
   graph(style=overlappingbar,noticks,klabel=keylab,key=upleft, $
   header='Histogram of Standardized Residuals') 2
   # nval 1 ng ; # yval 1 ng
   graph(max=1,min=-1,number=1,header='Correlogram of residuals', $
   hlabel='Lag',style=bar,shading=shad) 1
   # ser2 2 corrend+1
  spgraph(done)
 end do %i
 branch resids
 next
*
:corrs
*
 comp count = 0
 do %i=1,%DIMZ0
  do %j=1,%DIMZ0
   comp count = count + 1
   set ser1 / = UHAT(%i,t-%start-%LAG+1)
   set ser2 / = UHAT(%j,t-%start-%LAG+1)
   cross(noprint) ser1 ser2 / 0 corrend ser3
   set crosscorrel(count) 1 corrend+1 = ser3
  end do %j
 end do %i
 comp count = 0
 @%GRAFON%
 disp(store=header) 'Lags 1 to' %string(corrend)
 spgraph(vfields=%DIMZ0,hfileds=%DIMZ0,$
 header='Cross- and autocorrelograms of the residuals',Hlabel=header,$
 xlabels=alphalabel,ylabels=alphalabel)
  do %i=1,%DIMZ0
   do %j=1,%DIMZ0
    comp count = count+1
    set shad 1 corrend+1 = ABS(crosscorrel(count)) >=  up(1)
     graph(max=1,min=-1, $
     number=1,style=bar,noticks,scale=none,shading=shad) 1 %i %j
     # crosscorrel(count) 2 corrend+1
   end do %j
  end do %i
 spgraph(done)
 next
*
:cointvect
*
 dim graphlabel(%R)
 do %i=1,%R
  compute graphlabel(%i) = 'beta_' + %string(%i)
 end do %i
 select(prompt='Select one or more series',strings=graphlabel,status=%stop) $
 number selectlist
 if %stop == 0 .or. number==0 ; next
 do %i=1,number
  compute %whatg = selectlist(%i)
  set ser1 / = CI(%whatg,t-%start-%LAG+1)
  set ser2 / = BetaRk(%whatg,t-%start-%LAG+1)
  disp(store=header) 'beta' @-1 # %whatg @-1 '`  * Zk(t)'
  disp(store=header1) 'beta' @-1 # %whatg @-1 '`  * Rk(t)'
  @%GRAFON%
  spgraph(vfields=2)
   graph(header=header) 1
   # ser1
   graph(header=header1,samesize) 1
   # ser2 / 1
  spgraph(done)
 end do %i
 branch cointvect
 next
*
:vvectors
*
 dim graphlabel(%DimZ0)
 do %i=1,%DimZ0
  compute graphlabel(%i) = 'V_' + %string(%i)
 end do %i
 select(prompt='Select one or more series',strings=graphlabel,status=%stop) $
 number selectlist
 if %stop == 0 .or. number==0 ; next
 do %i=1,number
  compute %whatg = selectlist(%i)
  set ser1 / = ciu(%whatg,t-%start-%LAG+1)
  set ser2 / = BetaRku(%whatg,t-%start-%LAG+1)
  disp(store=header) 'V' @-1 # %whatg @-1 '`  * Zk(t)'
  disp(store=header1) 'V' @-1 # %whatg @-1 '`  * Rk(t)'
  @%GRAFON%
  spgraph(vfields=2)
   graph(header=header) 1
   # ser1
   graph(header=header1,samesize) 1
   # ser2 / 1
 spgraph(done)
 end do %i
 branch vvectors
 next
*
:savedat
 open(append) savedata catsdata.dat
 disp(unit=savedata) '/*'
 disp(unit=savedata) 'The model'
 disp(unit=savedata) 'beta (transposed)'
 write(unit=savedata,noskip,format='(2X8A8)') %betalab
 write(unit=savedata,noskip,format='(8F8.3)') %showbeta
 disp(unit=savedata) 'alpha'
 write(unit=savedata,noskip,format='(8F8.3)') %alpha1
 disp(unit=savedata) 'Sample' %datelabel(%start+%lag) '-' %datelabel(%end)
 disp(unit=savedata) '*/'
 loop
  disp(store=%prompt) 'Press 1 to save the ecm`s, 2 to save the residuals,'
  query(status=%stop,prompt=%prompt) %whatg
  if %stop == 0 ; break
  if %whatg == 1 ; begin
   disp(unit=savedata) '/*'
   disp(unit=savedata) 'The' %string(%R) 'ECM`s orderd by variable'
   disp(unit=savedata) '*/'
   write(unit=savedata,format='(8F9.5)') ci
   *** Added to save the Beta*Rk's as well:
   disp(unit=savedata) '/*'
   disp(unit=savedata) 'The' %string(%R) 'Beta*Rk`s orderd by variable'
   disp(unit=savedata) '*/'
   write(unit=savedata,format='(8F9.5)') betaRk


   cntrl(pause) 'The data has been saved in the file "catsdata.dat."'
  end
  else if %whatg == 2 ; begin
   disp(unit=savedata) '/*'
   disp(unit=savedata) 'The residuals orderd by variable'
   disp(unit=savedata) '*/'
   write(unit=savedata,format='(8F9.5)') uhat
   cntrl(pause) 'The data has been saved in the file "catsdata.dat."'
  end
  else ; next
 end loop
 close savedata
 next
*
:saveplot
 comp %gflag = %gflag+1
 next
*
:exit
/*
  Release memory
*/
 release ci ciu betark betarku fitmat uhat uhatmean Vhelp header
 release catserlabel alphalabel
*
 close plot
 return
*
end loop
*
END %GRAF%
*
PROCEDURE %RESCHK%
/*
Performs an analysis of the residuals in the model for any pi-matrix
*/
loc series ser1
loc sym  Suu rc covz
loc rect uhat uhat0 help1 help2 help3 help4
loc rect stat1 stat2 LM  param eigenvec
loc vect tLM seuhat skewness kurtosis eigenval Vhelp fr
loc vector[label] statlab1 statlab2
loc vect[series] endog
loc int corrend lag1 f1 f2
loc real N57 wilks
*
dim endog(%DIMZ0) rc(%DIMZ0,%DIMZ0) seuhat(%DIMZ0) tLM(10) fr(5)
dim stat1(%DIMZ0,6) stat2(%DIMZ0,3) Vhelp(%N)
dim statlab1(6) statlab2(3) skewness(%DIMZ0) kurtosis(%DIMZ0) param(2,4)
comp statlab1(1) =  'MEAN'
comp statlab1(2) =  'STD.DEV'
comp statlab1(3) =  'SKEWNESS'
comp statlab1(4) =  'KURTOSIS'
comp statlab1(5) =  'MAXIMUM'
comp statlab1(6) =  'MINIMUM'
comp statlab2(1) = 'ARCH(' + %LAG + ')'
comp statlab2(2) = 'Normality'
comp statlab2(3) = 'R-squared'
comp stat1    = %const(0)
comp stat2    = %const(0)
comp skewness = %const(0)
comp kurtosis = %const(0)
* comp  corrend = %N/4
if %n>20000
 compute corrend = 100
else
 compute corrend = %N/4
/*
 The covariance matrix for the endogenous series
*/
if %sflag == 1 ; begin
 smpl %start+%LAG %end
 do i = 1,%DIMZ0
  set endog(i) / = %Z0MAT(i,t-%start-%LAG+1)
 end do i
 *** Dec. 2000 fix here:
 compute %%sigma = %sigma
 VCV(centered,MATRIX=covz,noprint) / ; # endog
 compute %sigma = %%sigma

/*
 The parameters for the small sample correction in the test for normality
*/
 comp N57 = (%N1+5)*(%N1+7)
 comp param(1,1) = 3.0*(%N1**2+27*%N1-70)*(%N1+1)*(%N1+3)/((%N1-2)*N57*(%N1+9))
 comp param(1,2) = sqrt(2*(param(1,1)-1)) - 1
 comp param(1,3) = 1.0/sqrt(log(sqrt(param(1,2))))
 comp param(1,4) = sqrt((param(1,2)-1)*(%N1+1)*(%N1+3)/(12.0*(%N1-2)))
 comp param(2,1) = (%N1**2 + 15*%N1 - 4)*(%N1 - 3)*(%N1 + 1)
 comp param(2,2) = (%N1**2 + 27*%N1 -70)*(%N1 - 2)*N57/(6.0*param(2,1))
 comp param(2,3) = (%N1**2 +  2*%N1 - 5)*(%N1-7)*N57/(6.0*param(2,1))
 comp param(2,4) = (%N1**3 + 37*%N1**2 + 11*%N1 - 313)*N57/(12.0*param(2,1))
*
 comp %sflag=2
end
/*
 Calculate the residuals, covariance and correlation
*/
comp uhat = %Z0MAT - ((tr(%M10) - %pimat*tr(%M1k))*%iM11)*%Z1MAT - %pimat*%ZkMAT
ewise seuhat(i) = sqrt(%sigma(i,i))
ewise rc(i,j) = %sigma(i,j)/(seuhat(i)*seuhat(j))
/*
 Multivariate statistics
*/
comp tLM(1) = log(%det(%sigma))
comp tLM(2) = 1 - %dot(%sigma,inv(covz))/%DIMZ0
comp tLM(3) = %iN*(%DIMZ1*%DIMZ0+(%DIMZ0+%DIMZk-%R)*%R)*log(%N1) + tLM(1)
comp tLM(4) = %iN*(%DIMZ1*%DIMZ0+(%DIMZ0+%DIMZk-%R)*%R)*2*log(log(%N1)) + tLM(1)
/*
 Ljung-Box test
*/
comp tLM(5) = 0
comp help3 = inv(%sigma)
do %i = 2,corrend+1
 overlay uhat(1,1) with help1(%DIMZ0,%N-%i+1)
 overlay uhat(1,%i) with uhat0(%DIMZ0,%N-%i+1)
 comp help2 = (uhat0*tr(help1))*%iN
 comp LM = tr(help2)*help3*help2*help3
 comp tLM(5) = %trace(LM)*(1.0/(%N1-%i+1)) + tLM(5)
end do %i
comp tLM(5) = %N1*(%N1+2)*tLM(5)
/*
 Lagrange multiplier test for 1st and 4th order autocorrelation
*/
/*
 Parameters for small sample correction
*/
 comp fr(1) = %dimZ0**2
 comp fr(3) = %N1 - (%dimZ1+%R+%dimZ0) - 0.5
 comp fr(4) = sqrt((2*fr(1)-5.0)/(fr(1)**2 - 4))
 comp fr(5) = 0.5*fr(1)-1
 comp fr(2) = fr(3)/fr(4) - fr(5)
 comp f1 = fix(fr(1))
 comp f2 = fix(fr(2))
comp lag1 = 2
do %i=1,2
 comp help3 = tr(%beta1)*%ZkMAT
 overlay uhat(1,1) with help2(%DIMZ0,%N-lag1+1)
 dim help1(%DIMZ0,lag1-1)
 comp help1 = %const(0)
 @%CONCAT% help1 help2
 if %DIMZ1 > 0 ; @%CONCAT%(rows) help1 %Z1MAT
 @%CONCAT%(rows) help1 help3
 comp help2 = (help1*tr(help1))*%iN
 comp help3 = (%Z0MAT*tr(help1))*%iN
 comp Suu = %M00 - help3*inv(help2)*tr(help3)
 comp wilks = %det(Suu)/%det(%sigma)
 comp tLM(%i+5) = -1*fr(3)*log(wilks)
 comp tLM(%i+7) = fr(2)*(1.0-wilks**fr(4))/(fr(1)*wilks**fr(4))
 comp lag1 = 5
end do %i
/*
 Shenton-Bowman & "Doornik-Hansen" test for normality
*/
do %i=1,%DIMZ0
 ewise Vhelp(i) = uhat(%i,i)
 comp stat1(%i,1) = %sum(Vhelp)*%iN
end do %i
dim help1(%DIMZ0,%N) help2(%DIMZ0,1) help3(%DIMZ0,1) help4(%DIMZ0,1)
ewise help1(i,j) = (uhat(i,j)-stat1(i,1))/seuhat(i)
eigen rc eigenval eigenvec
ewise eigenval(i) = 1.0/(sqrt(eigenval(i)))
comp uhat0 = eigenvec*%diag(eigenval)*tr(eigenvec)*help1
do %i=1,%DIMZ0
 do %j=1,%N
  comp skewness(%i) = uhat0(%i,%j)**3 + skewness(%i)
  comp kurtosis(%i) = uhat0(%i,%j)**4 + kurtosis(%i)
  comp stat1(%i,3)  = help1(%i,%j)**3 + stat1(%i,3)
  comp stat1(%i,4)  = help1(%i,%j)**4 + stat1(%i,4)
 end do %j
end do %i
comp skewness = %iN*skewness
comp kurtosis = %iN*kurtosis
do %i=1,%DIMZ0
 comp stat1(%i,3)  = %iN*stat1(%i,3)
 comp stat1(%i,4)  = %iN*stat1(%i,4)
end do %i
*
ewise help2(i,j)  = param(1,4)*skewness(i)
ewise help3(i,j)  = param(2,2) + param(2,3)*(skewness(i)**2)
ewise help4(i,j)  = (kurtosis(i) -(skewness(i)**2) - 1)*param(2,4)
ewise skewness(i) = (param(1,3)*log(help2(i,1) + sqrt(help2(i,1)**2 + 1)))**2
ewise kurtosis(i) = $
((help4(i,1)/help3(i,1))**(1.0/3.0) -1 +1.0/(9*help3(i,1)))**2 *(9*help3(i,1))
comp tLM(10) = %sum(skewness) + %sum(kurtosis)
*
do %i=1,%DIMZ0
 comp help2(%i,1)  = param(1,4)*stat1(%i,3)
 comp help3(%i,1)  = param(2,2) + param(2,3)*stat1(%i,3)**2
 comp help4(%i,1)  = (stat1(%i,4) - stat1(%i,3)**2 - 1)*param(2,4)
end do %i
ewise skewness(i) = (param(1,3)*log(help2(i,1) + sqrt(help2(i,1)**2 + 1)))**2
ewise kurtosis(i) = $
((help4(i,1)/help3(i,1))**(1.0/3.0) -1 +1.0/(9*help3(i,1)))**2 *(9*help3(i,1))
do %i = 1,%DIMZ0
 comp stat2(%i,2) = skewness(%i) + kurtosis(%i)
end do %i
/*
 Write the Covariance matrix
*/
cntrl(clear)
disp(unit=&%copy) ' '
disp(unit=&%copy) 'RESIDUAL ANALYSIS '
disp(unit=&%copy) ' '
comp %outlab = 'output'
do %write = 1,%iwrite
 disp(unit=&%outlab) ' '
 disp(unit=&%outlab) 'Correlation matrix'
 write(unit=&%outlab,noskip,format='(8A10)') %alphalab
 write(unit=&%outlab,noskip,format='(8F10.6)') rc
 disp(unit=&%outlab) ' '
 disp(unit=&%outlab) 'Standard deviations of residuals'
 write(unit=&%outlab,format='(8F10.6)') seuhat
 comp %outlab = 'copy'
end do %write
disp(unit=&%copy) ' '
/*
 Write the multivariate statistics
*/
if %runtype==1 ; cntrl(pause,clear)
comp %outlab = 'output'
do %write = 1,%iwrite
 disp(unit=&%outlab) 'MULTIVARIATE STATISTICS'
 disp(unit=&%outlab) ' '
 disp(unit=&%outlab) 'LOG(DET(SIGMA))          =' tLM(1)
 disp(unit=&%outlab) 'INFORMATION CRITERIA: SC =' tLM(3)
 disp(unit=&%outlab) '                      HQ =' tLM(4)
 disp(unit=&%outlab) 'TRACE CORRELATION        =' tLM(2)
 disp(unit=&%outlab) ' '
 disp(unit=&%outlab) 'TEST FOR AUTOCORRELATION'
 disp(unit=&%outlab) 'L-B(' @-1 %string(corrend) @-1 '), CHISQ(' $
 @-1 ((%DIMZ0**2)*(corrend-%LAG+1)-%dimZ0*%R) @-1 ')' @21 '=' $
  ###.### tLM(5) @-1 ', p-val =' $
 #.## %chisqr(tLM(5),((%DIMZ0**2)*(corrend-%LAG+1)-%dimZ0*%R))
 disp(unit=&%outlab) 'LM(1),   CHISQ(' @-1 %string(f1) @-1 ')' @21 '=' $
 ###.### tLM(6) @-1 ', p-val =' #.## %chisqr(tLM(6),f1)
 disp(unit=&%outlab) 'LM(4),   CHISQ(' @-1 %string(f1) @-1 ')' @21 '=' $
 ###.### tLM(7) @-1 ', p-val =' #.## %chisqr(tLM(7),f1)
 disp(unit=&%outlab) ' '
 disp(unit=&%outlab) 'TEST FOR NORMALITY'
 disp(unit=&%outlab) 'CHISQ(' @-1 %string(2*%DIMZ0) @-1 ')' @21 '=' $
 ###.### tLM(10)  @-1 ', p-val =' #.## %chisqr(tLM(10),2*%DIMZ0)
 disp(unit=&%outlab) ' '
 comp %outlab = 'copy'
end do %write
/*
 Univariate residual analysis
*/
smpl %start+%LAG %end
do %i=1,%DIMZ0
 comp stat1(%i,2) = seuhat(%i)
 comp stat2(%i,3) = 1 - (stat1(%i,2)**2/covz(%i,%i))
 *
 set ser1 / = uhat(%i,t-%start-%LAG+1)
 *
 extremum(noprint) ser1
 comp stat1(%i,5)= %maximum
 comp stat1(%i,6)= %minimum
 *
 set ser1 / = uhat(%i,t-%start-%LAG+1)**2
 linreg(noprint) ser1 %start+2*%LAG %end
 # constant ser1{1 to %LAG}
 comp stat2(%i,1) = %nobs*%rsquared
end do %i
/*
 Write the univariate statistics
*/
if %runtype==1 ; cntrl(pause,clear)
comp %outlab = 'output'
do %write = 1,%iwrite
 disp(unit=&%outlab) 'UNIVARIATE STATISTICS'
 disp(unit=&%outlab) ' '
 write(unit=&%outlab,format='(6A12)',noskip) statlab1
 write(unit=&%outlab,noskip,format='(6F12.6)') stat1
 disp(unit=&%outlab) ' '
 write(unit=&%outlab,format='(5A12)',noskip) statlab2
 write(unit=&%outlab,format='(5F12.3)') stat2
 comp %outlab = 'copy'
end do %write
comp %outlab = 'output'
/*
 Release memory
*/
release Suu rc uhat uhat0 help1 help2 help3 stat1 stat2 LM
release tLM statlab1 statlab2 endog skewness kurtosis corrend
*
END %RESCHK%
*
PROCEDURE %COMPAN%
/*
This procedure calculates the eigenvalues
of the companion matrix for any pi-matrix
*/
loc rect VECM VARM BIGMAT ONES ZIPS IDEN gam gam1 roots
loc vect[complex] CROOTS
loc series xval yval xbord ybord
loc series[complex] root mod per
/*
 Calculate the coefficient matrix for the VAR-model
*/
cntrl(clear)
overlay %pimat(1,1) with VECM(%DIMZ0,%DIMZ0)
comp ONES = %identity(%LAG)
if %LAG > 1 ; begin
 comp gam = (tr(%M10) - %pimat*tr(%M1k))*%iM11
 overlay gam(1,1) with gam1(%DIMZ0,%DIMZ0*(%LAG-1))
 @%CONCAT% VECM gam1
 do %i=2,%LAG
  comp ONES(%i,%i) = -1
  comp ONES(%i,%i-1) = 1
 end do %i
 dim ZIPS(%DIMZ0,%DIMZ0*(%LAG-1))
 comp ZIPS = %const(0)
end
comp IDEN = %identity(%DIMZ0)
comp BIGMAT = %kroneker(ONES,IDEN)
if %LAG > 1 ; @%CONCAT% IDEN ZIPS
comp VARM = VECM*BIGMAT + IDEN
/*
 calculate the VAR(1)-rep (re-using the auxiliary matrices)
*/
if %LAG > 1 ; begin
 comp ONES = %identity(%LAG-1)
 comp IDEN = %identity(%DIMZ0)
 comp BIGMAT = %kroneker(ONES,IDEN)
 dim ZIPS(%DIMZ0*(%LAG-1),%DIMZ0)
 comp ZIPS = %const(0)
 @%CONCAT% BIGMAT ZIPS
 @%CONCAT%(ROWS) VARM BIGMAT
end
/*
 Calculate the eigenvalues of the companion-matrix
*/
eigen(cvalue=CROOTS,sort=absval) VARM
cset root 1 %DIMZ0*%LAG = croots(t)
polar root 1 %DIMZ0*%LAG mod per
dim roots(%DIMZ0*%LAG,4)
do %i=1,%DIMZ0*%LAG
 comp roots(%i,1) = %real(croots(%i))
 comp roots(%i,2) = %imag(croots(%i))
 comp roots(%i,3) = %real(mod(%i))
 comp roots(%i,4) = %real(per(%i))
end do %i
set xval 1 %DIMZ0*%LAG = roots(t,1)
set yval 1 %DIMZ0*%LAG = roots(t,2)
/*
 Create a unit circle to use in a scatter plot
*/
set xbord 1 100 = %real(%unit2(t,100))
set ybord 1 100 = %imag(%unit2(t,100))
*
scatter(style=symbols,extend=both, $
        header='The eigenvalues of the companion matrix') 2
# xbord ybord 1 100
# xval yval 1 %DIMZ0*%LAG
/*
 writing
*/
comp %outlab = 'output'
do %write=1,%iwrite
 disp(unit=&%outlab)  'The eigenvalues of the companion matrix'
 disp(unit=&%outlab) '      real   complex   modulus    argument'
 write(unit=&%outlab,format='(8F10.4)')  roots
 comp %outlab = 'copy'
end do %write
/*
 Release memory
*/
release VECM VARM BIGMAT ONES ZIPS IDEN gam gam1
*
END %COMPAN%
*
PROCEDURE %SHORTRUN%
/*
 Calculate the gamma matrix and t-values
*/
dec rect showgamma showtgamma detvarmat shortrunvarmat
*loc rect gamma tgamma showgamma showtgamma
loc rect gamma tgamma
loc rect  tgamma
loc sym Vgamma var2
loc int count
*
if %DIMZ1==0 ; begin
 messagebox 'YOU HAVE NO SHORT-RUN PARAMETERS'
 return
end
*
dim tgamma(%DIMZ0,%DIMZ1)
comp gamma = (tr(%M10) - %pimat*tr(%M1k))*%iM11
compute shortrunvarmat = gamma
comp var2 = %beta1*inv(%MQFORM(%Mkk,%beta1))*tr(%beta1)
comp Vgamma = %M11 - %M1k*var2*tr(%M1k)
comp Vgamma = inv(Vgamma)*%iN
do %i=1,%DIMZ0
 do %j=1,%DIMZ1
  comp tgamma(%i,%j) =  gamma(%i,%j)/sqrt(%sigma(%i,%i)*Vgamma(%j,%j))
 end do %j
end do %i
*
cntrl(clear)
disp(unit=&%copy) '  The short-run matrices'
disp(unit=&%copy) ' '
/*
 The lagged endogenous variables
*/
if %LAG > 1 ; begin
 disp(unit=&%copy) 'The lagged endogenous variables'
 dim %showlabel(%DIMZ0)
 comp count = 0
 do %i = 1,%LAG-1
  overlay gamma(1,count*%DIMZ0+1) with showgamma(%DIMZ0,%DIMZ0)
  overlay tgamma(1,count*%DIMZ0+1) with showtgamma(%DIMZ0,%DIMZ0)
  comp count = count+1
  comp %outlab ='output'
  do %write=1,%iwrite
   disp(unit=&%outlab) ' '
   disp(unit=&%outlab) 'Time: t-' @-1 %string(%i)
   disp(unit=&%outlab) ' '
   write(unit=&%outlab,format='(8A10)',noskip) %alphalab
   write(unit=&%outlab,format='(8F10.3)',noskip) showgamma
   disp(unit=&%outlab) ' '
   disp(unit=&%outlab) 't-values'
   write(unit=&%outlab,format='(8F10.3)',noskip) showtgamma
   disp(unit=&%outlab) ' '
   comp %outlab = 'copy'
  end do %write
  if %runtype==1 .and. %i < %LAG-1 ; cntrl(pause,clear)
 end do %i
end
/*
 The exogenous I(1) variables
*/
if %DIMX > 0 ; begin
 if %runtype==1 ; cntrl(pause,clear)
 disp(unit=&%copy) 'The differences of the exogenous I(1) variables'
 dim %showlabel(%DIMX)
 ewise %showlabel(j) = 'D' + %CATSERLAB(%DIMZ0+j)
 comp count = 0
 do %i = 0,%LAG-1
  overlay gamma(1,(%LAG-1)*%DIMZ0+count*%DIMX+1) with showgamma(%DIMZ0,%DIMX)
  overlay tgamma(1,(%LAG-1)*%DIMZ0+count*%DIMX+1) with showtgamma(%DIMZ0,%DIMX)
  comp count = count+1
  comp %outlab ='output'
  do %write=1,%iwrite
   disp(unit=&%outlab) ' '
   disp(unit=&%outlab) 'Time: t-' @-1 %string(%i)
   disp(unit=&%outlab) ' '
   write(unit=&%outlab,format='(8A10)',noskip) %showlabel
   write(unit=&%outlab,format='(8F10.3)',noskip) showgamma
   disp(unit=&%outlab) ' '
   disp(unit=&%outlab) 't-values'
   write(unit=&%outlab,format='(8F10.3)',noskip) showtgamma
   disp(unit=&%outlab) ' '
   comp %outlab = 'copy'
  end do %write
  if %runtype==1 .and. %i < %LAG-1 ; cntrl(pause,clear)
 end do %i
end
/*
 The exogenous I(0) variables
*/
if %DIMD > 0 ; begin
 if %runtype==1 ; cntrl(pause,clear)
 disp(unit=&%copy) 'The exogenous I(0) variables'
 dim %showlabel(%DIMD)
 ewise %showlabel(j) =  %CATSERLAB(%DIMZ0+%DIMX+j)
 overlay gamma(1,(%LAG-1)*%DIMZ0+%LAG*%DIMX+1) with showgamma(%DIMZ0,%DIMD)
 overlay tgamma(1,(%LAG-1)*%DIMZ0+%LAG*%DIMX+1) with showtgamma(%DIMZ0,%DIMD)
 comp %outlab ='output'
 do %write=1,%iwrite
  disp(unit=&%outlab) ' '
  write(unit=&%outlab,format='(8A9)',noskip) %showlabel
  write(unit=&%outlab,format='(8F9.3)',noskip) showgamma
  disp(unit=&%outlab) ' '
  disp(unit=&%outlab) 't-values'
  write(unit=&%outlab,format='(8F9.3)',noskip) showtgamma
  disp(unit=&%outlab) ' '
  comp %outlab = 'copy'
 end do %write
end
/*
* The deterministic variables *
*/
if %SEA > 0 .OR. %DETER > 2 ; begin
 if %runtype==1 ; cntrl(pause,clear)
 disp(unit=&%copy) 'The deterministic variables'
 comp count = FIX(%IF(%SEA>0,%SEA,0)) + FIX(%IF(%DETER>2,1,0))
 overlay gamma(1,(%LAG-1)*%DIMZ0+%LAG*%DIMX+%DIMD+1) $
  with showgamma(%DIMZ0,count)
 overlay tgamma(1,(%LAG-1)*%DIMZ0+%LAG*%DIMX+%DIMD+1) $
  with showtgamma(%DIMZ0,count)
 dim %showlabel(count)
 do %i=1,%SEA
  comp %showlabel(%i) =  'SEA(' + %string(%i) +')'
 end do %i
 if %DETER > 2 ; comp %showlabel(count) = 'CONST'
 comp %outlab ='output'
 do %write=1,%iwrite
  disp(unit=&%outlab) ' '
  write(unit=&%outlab,format='(12A7)',noskip) %showlabel
  write(unit=&%outlab,format='(12F7.3)',noskip) showgamma
  compute detvarmat = showgamma
  disp(unit=&%outlab) ' '
  disp(unit=&%outlab) 't-values'
  write(unit=&%outlab,format='(12F7.3)',noskip) showtgamma
  disp(unit=&%outlab) ' '
  comp %outlab = 'copy'
 end do %write
end
/*
 Release memory
*/
*release gamma tgamma Vgamma showgamma showtgamma var2 count
release gamma tgamma Vgamma  var2 count
*
END %SHORTRUN%
*
PROCEDURE %MAREP%
/*
 Calculate the moving average parameters beta-,alpha-orthogonal
 the impact matrix C(1) and the deterministic trend components.
*/
loc rect gam1 Gamma C_1 tC_1 mu tau help1 Q Q1
loc rect alfao betao betaX betat ksi load ones vmat variance1
loc sym CSC Mom
loc int no count
*
cntrl(clear)
*
dim betaX(%DIMZ0,%R) betat(1,%R) mu(%DIMZ0,1) %showlabel(%DIMZ0)
dim tC_1(%DimZ0,%DimZ0)
comp mu = %const(0)
ewise %showlabel(%i) = %betalab(%i)
ewise betaX(i,j) = %beta1(i,j)
if %DETER == 2 .OR. %DETER == 4
do %j=1,%R
 comp betat(1,%j) = %beta1(%DIMZk,%j)
end do %j
*
if %R == %DimZ0 ; begin
 comp C_1 = %scalar(-1.0)*inv(%alpha1*tr(betaX))
 disp 'Note: The long-run matrices have full rank'
 disp ' '
 branch fullrank
end
*
@%ORTHO% %alpha1 alfao
@%ORTHO% betaX betao
/*
 NOTE: NEW DEFINITION OF ALPHAO
*/
comp alfao = tr(alfao)
/*
 Calculating the impact matrix for the MA-rep
*/
comp gam1 = (tr(%M10) - %pimat*tr(%M1k))*%iM11
comp Gamma = %Identity(%DimZ0)
if %LAG > 1 ; begin
 comp count = 0
 dim help1(%DimZ0,%DimZ0)
 do %i=1,%LAG-1
  ewise help1(i,j) = gam1(i,count*%DimZ0+j)
  comp Gamma = Gamma - help1
  comp count = count + 1
 end do %i
end
*
comp C_1 = betao*inv(alfao*Gamma*betao)*alfao
comp CSC = C_1*%Sigma*tr(C_1)
/*
Calculating "t-values" for C_1
*/
comp vmat = tr(%beta1)*%ZkMAT
comp Q = (tr(C_1)*tr(Gamma)-%Identity(%DimZ0))*%alpha1*inv(tr(%alpha1)*%alpha1)
*
if %LAG > 1 ; begin
 dim Ones(1,%LAG-1) help1(%DimZ0*(%LAG-1),%N)
 comp Ones = %const(1)
 comp Q1 = tr(C_1)*%Kroneker(Ones,%Identity(%DimZ0))
 @%CONCAT% Q Q1
 ewise help1(i,j) = %Z1MAT(i,j)
 @%CONCAT%(ROWS) vmat help1
end
*
comp Mom = vmat*tr(vmat)
*
if %DETER>2 .or. %SEA>0 .or. %DimD>0 .or. %DimX>0 ; begin
 dim help1(%ROWS(%Z1MAT)-%DimZ0*(%LAG-1),%N)
 ewise help1(i,j) = %Z1MAT(%DimZ0*(%LAG-1)+i,j)
 comp Mom = Mom - vmat*tr(help1)*inv(help1*tr(help1))*help1*tr(vmat)
end if
*
compute variance1 = %MQFORMDIAG(inv(Mom),tr(Q))
*
do %i=1,%DimZ0
 do %j=1,%DimZ0
  comp tC_1(%i,%j) = C_1(%i,%j)/sqrt(CSC(%i,%i)*variance1(%j,1))
 end do %j
end do %i
/*
 Calculating the trend coefficients
*/
if %DETER > 2
do %i=1,%DIMZ0
  comp mu(%i,1) = gam1(%i,%DIMZ1)
end do %i
comp tau = C_1*mu
if %DETER == 4 ; begin
 comp help1 = C_1*Gamma - %Identity(%DIMZ0)
 comp tau = tau + help1*betaX*inv(tr(betaX)*betaX)*tr(betat)
end
comp tau = tr(tau)
/*
  Normalization of alfao
*/
comp ksi = %identity(%DIMZ0-%R)
disp 'Re-normalisation of the orthogonal complement to alpha'
disp 'The orthogonal complement to alpha (transposed)'
write(noskip,format='(8A10)') %showlabel
write(format='(8F9.4)') alfao
do %i=1,%DIMZ0-%R
 :error
 disp(store=%prompt) $
 'Input the element of vector' %string(%i) 'to use for transformation'
 query(status=%stop,prompt=%prompt) no
 if %stop == 0 .or. no == 0 ; break
 if no .lt. 1 .or. no .gt. %DIMZ0 ;  begin
  cntrl(pause) 'Sorry, wrong input. Please re-enter'
  branch error
 end
 if alfao(%i,no) == 0 ; begin
  comp %prompt = $
  'Sorry, normalizing on 0 in row_' + %string(%i) + '. Please re-enter.'
  messagebox &%prompt
  branch error
 end
 comp ksi(%i,%i) = alfao(%i,no)
end do %i
*
comp alfao = inv(ksi)*alfao
comp load = betao*inv(alfao*Gamma*betao)
*
cntrl(clear)
comp %outlab = 'output'
do %write=1,%iwrite
 disp(unit=&%outlab) 'THE MA-REPRESENTATION AND DECOMPOSITION OF THE TREND'
 disp(unit=&%outlab) ' '
 disp(unit=&%outlab) 'The orthogonal complement to alpha'
 disp(unit=&%outlab) '(The coefficients to the common trends (transposed))'
 write(unit=&%outlab,noskip,format='(2X8A10)') %showlabel
 write(unit=&%outlab,noskip,format='(8F10.3)') alfao
 disp(unit=&%outlab) ' '
 disp(unit=&%outlab) 'The loadings to the common trends'
 write(unit=&%outlab,noskip,format='(8F10.3)') load
 disp(unit=&%outlab) ' '
 comp %outlab = 'copy'
end do %write
*
if %runtype==1 ; cntrl(pause,clear)
*
comp %outlab = 'output'
do %write=1,%iwrite
 disp(unit=&%outlab) 'The impact-matrix for the MA-representation: C(1)'
 write(unit=&%outlab,noskip,format='(2X8A10)') %showlabel
 write(unit=&%outlab,noskip,format='(8F10.3)') C_1
 disp(unit=&%outlab) ' '
 disp(unit=&%outlab) '"t-values" for C(1)'
 write(unit=&%outlab,noskip,format='(8F10.3)') tC_1
 disp(unit=&%outlab) ' '
 comp %outlab = 'copy'
end do %write
*
if %runtype==1 ; cntrl(pause,clear)
*
comp %outlab = 'output'
do %write=1,%iwrite
 disp(unit=&%outlab) "The long-run covariance matrix: C(1)*Sigma*C(1)'"
 write(unit=&%outlab,noskip,format='(2X8A10)') %showlabel
 write(unit=&%outlab,noskip,format='(8F10.4)') CSC
 disp(unit=&%outlab) ' '
 *
 if %DETER > 2 ; begin
  disp(unit=&%outlab) 'The linear trends in the levels'
  write(unit=&%outlab,noskip,format='(2X8A10)') %showlabel
  write(unit=&%outlab,noskip,format='(8F10.4)') tau
  disp(unit=&%outlab) ' '
 end
 comp %outlab = 'copy'
end do %write
*
branch exit
*
:fullrank
*
comp %outlab = 'output'
do %write=1,%iwrite
 disp(unit=&%outlab) 'The impact-matrix for the MA-representation: C(1)'
 write(unit=&%outlab,noskip,format='(2X8A10)') %showlabel
 write(unit=&%outlab,noskip,format='(8F10.3)') C_1
 disp(unit=&%outlab) ' '
 comp %outlab = 'copy'
end do %write
*
:exit
/*
 Release memory
*/
release gam1 Gamma C_1 tC_1 mu tau help1 Q Q1
release alfao betao betaX betat ksi load ones vmat variance1 CSC Mom
release no count
*
END %MAREP%

