Ignacio Diaz-Emparanza 1.0 2007-06-05 Local Linear Trend Model

if (irreg!=1&&irreg!=0) || (level!=1&&level!=0) || (slope!=1&&slope!=0)
  funcerr "irreg, level and slope must be 0 or 1"
end if
if $pd>1
  print "LLTestim warning: your data are seasonal, you should use the BSMestim function"
end if
matrix fijos = { irreg, level, slope }
matrix theta = { -0.5, -1.5}
matrix at
series V = 0
series F = 0
scalar logLc=0
matrix pstar
scalar sigma_e=1
matrix bT
set bfgs_toler 1.E-2
M = BFGSmax(theta, "LLT(&theta, y, sigma_e, &bT, &at, &pstar, &V, &F, &logLc, sigmatol, 1, fijos)")
scalar sstart = int(min(t))
scalar send = int(max(t))
scalar T=send-sstart+1
scalar sstar = (1/T)*(sum((V^2)/F)
scalar q1=exp(2*theta[1])*sstar
scalar q2=exp(2*theta[2])*sstar
matrix qp = { sstar, q1, q2 }
matrix conc=imaxr(qp)
scalar concent = conc[1]
set bfgs_toler default
matrix theta = { -0.5, -1.5 }
M = BFGSmax(theta, "LLT(&theta, y, sigma_e, &bT, &at, &pstar, &V, &F, &logLc, sigmatol, concent, fijos)")
scalar sstar = (1/T)*(sum((V^2)/F)
scalar q1=exp(2*theta[1])
scalar q2=exp(2*theta[2])
series LLTtrend = kf_smooth(pstar, &at, bT)
series LLTslope = at[2,]
list compo = LLTtrend LLTslope
printf "\nLocal Linear Trend Model estimation:\n"
printf "-----------------------------------------\n"
if concent=1
  printf "    sigma*=\t\t Var(eps)=%8.6E\n", sstar
  printf "    q1=%8.5f,\t Var(eta)=%8.6E\n", q1, q1*sstar
  printf "    q2=%8.5f,\t Var(xi)=%8.6E\n", q2, q2*sstar
  matrix q = {1, q1, q2}'
else
  if concent=2
    printf "    q1=%8.5f,\t Var(eps)=%8.6E\n", q1, q1*sstar
    printf "    sigma*=\t\t Var(eta)=%8.6E\n", sstar
    printf "    q2=%8.5f,\t Var(xi)=%8.6E\n", q2, q2*sstar
    matrix q = {q1, 1, q2}'
  else
    if concent=3
      printf "q1=%8.5f,\t Var(eps)=%8.6E\n", q1, q1*sstar
      printf "q2=%8.5f,\t Var(eta)=%8.6E\n", q2, q2*sstar
      printf "sigma*=\t\t Var(xi)=%8.6E\n", sstar
      matrix q = {q1, q2, 1}'
    end if
  end if
end if
printf "------------------------------------------\n \n"

/*
Measurement equation:
y(t) = Z[t,]*alpha(t)+e(t)	(1.1a)
State transition:
alpha(t)=bT*alpha(t-1)+w(t);  (1.2a)
bT is for "bold T" and w(t)=R(t)*eta(t) in Harvey 1990
a(t) is the estimator of alpha(t)
Parameters:
y       = observable series
a0      = m x 1 vector, prior a(0)
p0      = m x m matrix, prior p(0)=var(a(0))
bT      = m x m matrix (transition matrix)
Z       = T x m matrix
Sigma_w = m x m symmetric matrix of variance of w(t), fixed for all t
sigma_e = scalar variance of e(t), fixed for all t
at      = m x T  matrix (output) with the estimated states
pstar   = m^2 x T matrix (output)
*/
#
# Forward solution
#
scalar T = rows(Z)
scalar m = rows(a0)
matrix at_t = a0
matrix pt_t = p0
matrix at = zeros(m,T)
# printf "\n...Filtering...\n"
loop for i=1..T --quiet
  # Prediction equations
  # eq. (2.2a)
  matrix at_t = bT*at_t
  # eq. (2.2b)
  matrix pt_t1 = qform(bT,pt_t)+Sigma_w
  # eq. (2.3c)
  matrix zt = Z[$i,]
  matrix H = pt_t1*zt'
  matrix f = zt*H + sigma_e
  if i>1
    matrix pstar_t = (bT*pt_t)' inv(pt_t1)
  endif
  #Updating equations
  # eq. (2.4a)
  genr V[$i] = y[$i] - zt*at_t
  matrix at_t = at_t + H*(V[$i]/f)
  genr F[$i] = f
  # eq (2.4b)
  matrix pt_t = pt_t1 - H*H' * (1/f)
  matrix at[,$i]=at_t
  if i>1
    if i=2
      matrix pstar=vec(pstar_t)
    else
      matrix pstar=pstar~vec(pstar_t)
    endif
  endif
end loop
#Concentrated Log-likelihood fuction
scalar logLc = -(T/2)*(log(2*pi)+1)-(1/2)*sum(log(F))-(T/2)*log((1/T)*sum((V^2)/F)
series filtered = at[1,]
#  printf "\nFilter done\n"

/*
Fixed-interval smoothing
The matrix pt is not used here, but could be used if
one wants to calculate confidence intervals for the
at estimators.
*/
scalar m = rows(at)
scalar T = cols(at)
#printf "\n...Smoothing...\n"
scalar T1=T-1
loop for i=1..T1 --quiet
  scalar j=T-i
  matrix pstar_t = mshape(pstar[,j],m,m)
  # eq. (2.9a)
  matrix at[,j] += pstar_t*(at[,(j+1)]-bT*at[,j])
end loop
series ret = at[1,]
#printf "\nSmoothing done\n"

if concent>4
  funcerr "concent must be 1, 2, or 3"
end if
genr time
scalar sstart = int(min(time))
scalar send = int(max(time))
scalar T=send-sstart+1
matrix bT = { 1, 1; 0, 1 }
scalar m = cols(bT)
matrix Z = ones(T,1) ~ zeros(T,1)
matrix a0 = y[sstart] * ones(2,1)
matrix p0 = 400000*I(2)
matrix Sigma_w = zeros(2,2)
if concent=1
  scalar sigma_e=1
  scalar tmp = exp(2*param[1])*fixed[2]
  Sigma_w[1,1] = (tmp>sigmatol) ? tmp : 0
  matrix param[1] = (tmp>sigmatol) ? param[1] : -500
  scalar tmp = exp(2*param[2])*fixed[3]
  Sigma_w[2,2] = (tmp>sigmatol) ? tmp : 0
  matrix param[2] = (tmp>sigmatol) ? param[2] : -500
else
  if concent=2
    scalar tmp = exp(2*param[1])*fixed[1]
    scalar sigma_e=(tmp>sigmatol) ? tmp : 0
    matrix param[1] = (tmp>sigmatol) ? param[1] : -500
    Sigma_w[1,1] = 1
    scalar tmp = exp(2*param[2])*fixed[3]
    Sigma_w[2,2] = (tmp>sigmatol) ? tmp : 0
    matrix param[2] = (tmp>sigmatol) ? param[2] : -500
  else
    if concent=3
      scalar tmp = exp(2*param[1])*fixed[1]
      scalar sigma_e=(tmp>sigmatol) ? tmp : 0
      matrix param[1] = (tmp>sigmatol) ? param[1] : -500
      Sigma_w[2,2] = 1
      scalar tmp = exp(2*param[2])*fixed[2]
      Sigma_w[1,1] = (tmp>sigmatol) ? tmp : 0
      matrix param[2] = (tmp>sigmatol) ? param[2] : -500
    endif
  endif
endif
series filtered = kf_filt(y, a0, p0, bT, Z, Sigma_w, sigma_e, &at, &pstar, &V, &F, &logLc)