Estimating a TVP Model With the FKF Package in R

This example shows how to estimate a time-varying parameter model using the FKF package in R. That package is no longer available from CRAN because a message to the maintainer’s email bounced, but it’s a small package with stable code, so there is no obvious reason you shouldn’t continue to use it. In any event, much of the setup information will carry over to other packages.

The model we want to estimate

The model is \[y_{t} = \beta_{0} + \beta_{1} w_{t} + \beta_{2} v_{t} + \beta_{3t} x_{t} + \beta_{4t} z_{t} + \varepsilon_{yt},\] with \(\beta_{0}\), \(\beta_{1}\), and \(\beta_{2}\) constant over the sample, \(var\left( \varepsilon_{y} \right) = g^{2}\), \[\beta_{3t} = \beta_{3, t-1} + \eta_{3t},\] \[\beta_{4t} = \beta_{4, t-1} + \eta_{4t},\] \(var\left( \eta_{3} \right) = h_{3}\), and \(var\left( \eta_{4} \right) = h_{4}\).

State space form: state equation

Defining the state vector to be \[\alpha_{t} = \left[\begin{array}{c} \beta_{3t}\\ \beta_{4t}\\ \end{array}\right],\] the state equation is \[\left[\begin{array}{c} \beta_{3t}\\ \beta_{4t}\\ \end{array}\right] = \left[\begin{array}{c} \beta_{3, t-1}\\ \beta_{4, t-1}\\ \end{array}\right] + \left[\begin{array}{cc} h_{3} & 0\\ 0 & h_{4}\\ \end{array}\right] \left[\begin{array}{c} \eta_{3t}\\ \eta_{4t}\\ \end{array}\right].\]

State space form: measurement equation

Let’s rewrite the estimation equation in a way that isolates the state vector: \[y_{t} = \beta_{0} + \beta_{1} w_{t} + \beta_{2} v_{t} + \left[\begin{array}{cc} x_{t} & z_{t}\\ \end{array}\right] \left[\begin{array}{c} \beta_{3t}\\ \beta_{4t}\\ \end{array}\right] + g \varepsilon_{t}.\] Note that we’ve changed the error term such that \(var\left( \varepsilon \right) = 1\). The only reason we’ve done that is because it is necessary in order to fit it into the form required by FKF to do the estimation. \(\varepsilon_{t}\) and \(\varepsilon_{yt}\) (the error in the original equation) are different: \(\varepsilon_{yt} = g \varepsilon_{t}\).

FKF required form

We now need to put the state and measurement equations in the form required by FKF. The state equation has to be in the form \[\alpha_{t} = d_{t} + T_{t} \alpha_{t} + H_{t} \eta_{t},\] where the elements of \(\eta_{t}\) are all distributed \(N \left( 0,1 \right)\). We have that if \[d_{t} = 0\] \[T_{t} = I_{2}\] \[H_{t} = \left[\begin{array}{cc} h_{3} & 0\\ 0 & h_{4}\\ \end{array}\right].\] We do not need to specify anything about \(\alpha_{t}\), because the Kalman filter takes care of the state vector, or \(\eta_{t}\), because that is a normally distributed random error.

The measurement equation has to take the form \[y_{t} = c_{t} + Z_{t} \alpha_{t} + G_{t} \varepsilon_{t}.\] Define \[c_{t} = \beta_{0} + \beta_{1} y_{t-1} + \beta_{2} y_{t-2}\] \[Z_{t} = \left[\begin{array}{cc} x & z\\ \end{array}\right]\] \[G_{t} = g.\]

Estimation

The unknown parameters are \[\theta = \left[\begin{array}{c} \beta_{0}\\ \beta_{1}\\ \beta_{2}\\ h_{3}\\ h_{4}\\ g\\ \end{array}\right]\]

If we knew \(\theta\), the Kalman filter could be used to compute the state vector (\(\beta_{3}\) and \(\beta_{4}\)) at each \(t\). That facilitates estimation using the following algorithm:

1. Make an initial guess for \(\theta\).
2. Calculate \(c_{t}\), \(Z_{t}\), \(G_{t}\), \(T_{t}\), and \(H_{t}\) using \(\theta\).
3. Call the fkf function to evaluate the likelihood function.
4. Update \(\theta\).
5. Repeat steps 2-4 until you are satisfied that the current choice of \(\theta\) maximizes the likelihood.

Using the fkf function

We will demonstrate the use of the fkf function using generated data. There will be 500 observations, with large changes in \(\beta_{3}\) and \(\beta_{4}\) occurring at observations 150 and 350. We will see how well the TVP model is able to accommodate those changes.

Some notes on the notation in the fkf documentation:

• n: number of observations on y
• d: number of y variables (equal to one if you have a linear regression)
• m: number of state variables

Generate the data:

set.seed(100)
w <- rnorm(500)
v <- rnorm(500)
x <- rnorm(500)
z <- rnorm(500)
y <- local({
  b0 <- 0.2
  b1 <- 1.0
  b2 <- 1.0
  b3 <- c(rep(0.0,150), rep(0.5,200), rep(1.0,150))
  b4 <- b3
  g <- 1.0
  b0 + b1*w + b2*v + b3*x + b4*z + g*rnorm(500)
})

Now do the estimation. We need to write a function that takes a guess for \(\theta\) and returns the negative of the loglikelihood function.

objfnc <- function(theta) {
  # Destructure theta
  b0 <- theta[1]
  b1 <- theta[2]
  b2 <- theta[3]
  h3 <- theta[4]
  h4 <- theta[5]
  g <- theta[6]
  
  # Compute the matrices for the measurement equation
  ct <- matrix(b0 + b1*w + b2*v, nrow=1)
  zt <- array(dim=c(1,2,500))
  zt[1, 1,] <- x
  zt[1, 2,] <- z
  GGt <- array(dim=c(1,1,1))
  GGt[1,1,1] <- g
  
  # Compute the matrices for the state equation
  dt <- matrix(c(0,0), nrow=2)
  Tt <- array(dim=c(2,2,1))
  Tt[ , , 1] <- diag(2)
  HHt <- array(dim=c(2,2,1))
  HHt[1,1,1] <- h3
  HHt[2,2,1] <- h4
  HHt[1,2,1] <- HHt[2,1,1] <- 0

  # Initial values of the state vector
  a0 <- c(0,0)
  
  # Variances of a0
  # Follows the example in the package
  P0 <- matrix(1000000, nrow=2, ncol=2)
  
  # The data
  yt <- matrix(y, nrow=1)

  # Run the Kalman filter
  library(FKF)
  kf <- fkf(a0, P0, dt, ct, Tt, zt, HHt, GGt, yt)
  
  # Return the negative of the loglikelihood
  return(-kf$logLik)
}

And now maximize the likelihood function. We can use optim because objfnc returns the negative of the likelihood function.

fit <- optim(c(0,0,0,1,1,1), objfnc, control=list(maxit=10000)) 
print(fit$par)

We now have the estimated parameters, but we still need to recover the filtered state variables (\(b_{3t}\) and \(b_{4t}\)).

b0 <- fit$par[1]
b1 <- fit$par[2]
b2 <- fit$par[3]
h3 <- fit$par[4]
h4 <- fit$par[5]
g <-  fit$par[6]

ct <- matrix(b0 + b1*w + b2*v, nrow=1)
zt <- array(dim=c(1,2,500))
zt[1, 1,] <- x
zt[1, 2,] <- z
GGt <- array(dim=c(1,1,1))
GGt[1,1,1] <- g

dt <- matrix(c(0,0), nrow=2)
Tt <- array(dim=c(2,2,1))
Tt[ , , 1] <- diag(2)
HHt <- array(dim=c(2,2,1))
HHt[1,1,1] <- h3
HHt[2,2,1] <- h4
HHt[1,2,1] <- HHt[2,1,1] <- 0

a0 <- c(0,0)
P0 <- matrix(1000000, nrow=2, ncol=2)
yt <- matrix(y, nrow=1)
kf <- fkf(a0, P0, dt, ct, Tt, zt, HHt, GGt, yt)
pdf("tvp-coefficients.pdf")
  plot(kf$att[1,], type="l", main="Estimates of b3")
  plot(kf$att[2,], type="l", main="Estimates of b4")
dev.off()