*
* TVARYING.RPF
* RATS Version 8, User's Guide, Example from Section 7.14.2,7.14.3
* Typo in decayfac (wrong sign) corrected 10/2007. Thanks to
* Rodolfo Mendez for pointing this out.
*
cal(q) 1959
open data haversample.rat
data(format=rats) 1959:1 2006:4 gdph fm1

log gdph
log fm1

compute kstart=1960:1,kend=2006:4

equation kalmeq gdph
# constant gdph{1 to 4} fm1{1 to 4}
declare vector b(9) bmean(9)
ewise bmean(i)=(i==2)
associate(perm) kalmeq b
*
linreg gdph
# constant gdph{1 to 4}
compute gseesq=.9*%seesq
linreg fm1
# constant fm1{1 to 4}
compute mseesq=.9*%seesq
*
*  Set up the system. Use GSEESQ as the measurement equation variance
*
system kalmeq
kfset(constant,noscale,likelihood=likely) xxx
# gseesq
tvarying tvx
end(system)
*
* The hyperparameters (PIx) below are the following:
*   PI5 = Overall tightness. This corresponds to TIGHTNESS**2 on
*         SPECIFY.
*   PI2 = Relative tightness on other variables. This corresponds
*         to w**2 on the symmetric prior.
*   PI3 = Relative tightness on the constant. (No correspondence).
*   PI7 = Relative tightness on time variation.
*   PI8 = Shrinkage factor toward mean (between 0 and 1 , 1 = no shrinkage)
*
*  This uses a HARMONIC decay with DECAY=.5
*
compute pi5 = .20**2
compute pi2 = .5**2
compute pi3 = 10000.
compute pi7 = 0.00000001
compute pi8 = 1.00
*
dim xxx(9,9) tvx(9,9)
*
compute likely=%const(0.0)
compute b=bmean
compute xxx=%const(0.0)
compute xxx(1,1)=pi5*pi3*gseesq
do i=1,4
  compute decayfac=i**(-2*.5)
  compute xxx(i+1,i+1)=pi5*decayfac
  compute xxx(i+5,i+5)=pi5*pi2*decayfac*gseesq/mseesq
end do i
compute tvx=pi7*xxx
*
do time=kstart,kend
  if time==kstart
     kalman(start=time)
  else {
     compute b=pi8*b+(1-pi8)*bmean
     compute xxx=(pi8**2)*xxx
     kalman(print=(time==kend))
     }
end do time
*
*  Compute the concentrated log likelihood
*
compute nobs=kend-kstart+1
compute loglikely=-.5*(likely(1,1)+nobs*log(likely(2,1)/nobs))

compute pi5 = .20**2 , pi2 = .5**2 , pi3 = 10000.
compute pi7 = 0.00000001 , pi8 = 1.00
nonlin pi8 pi7
dec real loglikely
dim xxx(9,9) tvx(9,9)
find(method=genetic,iterations=200) maximum loglikely {
   if pi7<0.or.pi8>1.00 {
      compute loglikely=%na
      next
   }
   compute likely=%const(0.0)
   compute b=bmean
   compute xxx=%const(0.0)
   compute xxx(1,1)=pi5*pi3*gseesq
   do i=1,4
     compute decayfac=i**(-2*.5)
     compute xxx(i+1,i+1)=pi5*decayfac
     compute xxx(i+5,i+5)=pi5*pi2*decayfac*gseesq/mseesq
   end do i
   compute tvx=pi7*xxx
*
   do time=kstart,kend
      if time==kstart
	      kalman(start=time)
      else {
	       compute b=pi8*b+(1-pi8)*bmean
	       compute xxx=(pi8**2)*xxx
	       kalman
	   }
   end do time
   compute nobs=kend-kstart+1
   compute loglikely=-.5*(likely(1,1)+nobs*log(likely(2,1)/nobs))
}
end find