Jumpy TVPs

[fipelle]

Hi,

I have written the following code to estimate a univariate linear regression with time-varying parameters (TVP), defining the dynamics of the betas as integrated random walks. In theory, this should result in a relatively smooth outcome. However, the results are very jumpy. Am I setting up the function incorrectly? If not, how can I change the optimizer to something else, perhaps using an EM algorithm, Newton's method or decrease the tolerance?

tvp_regression <- function(y, x) {

    # Define `T`
    T_len <- length(y)

    # Define the 1x2xT `Z_array`
    Z_array <- array(0, dim = c(1, 2, T_len))
    Z_array[1, 1, ] <- x                      # Link beta_t with x_t
    # Z_array[1, 2, ] remains 0, ensuring slope_t does not affect observations
    
    # Define the 2x2 matrices for the smooth trend dynamics
    T_matrix <- matrix(c(1, 1, 0, 1), nrow=2, byrow=TRUE)
    Q_matrix <- matrix(c(0, 0, 0, NA), nrow=2, byrow=TRUE)

    # Define the initial state mean (a1) and covariance (P1) for diffuse states
    a1 <- c(0, 0)               # Initial mean
    P1 <- matrix(0, 2, 2)       # Initial state covariance (zero for diffuse initialization)
    P1inf <- diag(2)            # Diffuse indicator matrix (identity matrix for both states)

    # Define the state-space model
    model <- SSModel(
        y ~ -1 + SSMcustom(
                Z = Z_array,    # Observation matrix: y(t) = beta(t)*x(t) + noise
                T = T_matrix,   # State transition matrix
                R = diag(2),    # State noise selection matrix
                Q = Q_matrix,   # State noise covariance matrix
                a1 = a1,        # Initial mean
                P1 = P1,        # Initial covariance
                P1inf = P1inf   # Selection matrix for the exact diffuse initialisation 
            ),
            H = NA              # Observation noise variance to be estimated
    )
    
    # Estimate the parameters via MLE
    fit <- fitSSM(model, inits=c(0.1, 0.1), method="BFGS")
    
    # Get filtered estimates (the smoother is not needed here, as I only need the last value)
    filtered <- KFS(fit$model, smoothing="none")

    # Extract the most recent regression coefficient
    latest_beta <- filtered$att[T_len, 1]
    
    # Return output
    return(latest_beta)

[helske]

Using simulated data this seems to work as expected:

n <- 100
set.seed(1)
x <- 2 + rnorm(n)
beta <- cumsum(cumsum(rnorm(n, sd = 0.1)))

y <- beta * x + rnorm(n)
ts.plot(beta)
ts.plot(y)


# Define `T`
T_len <- length(y)

# Define the 1x2xT `Z_array`
Z_array <- array(0, dim = c(1, 2, T_len))
Z_array[1, 1, ] <- x                      # Link beta_t with x_t
# Z_array[1, 2, ] remains 0, ensuring slope_t does not affect observations

# Define the 2x2 matrices for the smooth trend dynamics
T_matrix <- matrix(c(1, 1, 0, 1), nrow=2, byrow=TRUE)
Q_matrix <- matrix(c(0, 0, 0, NA), nrow=2, byrow=TRUE)

# Define the initial state mean (a1) and covariance (P1) for diffuse states
a1 <- c(0, 0)               # Initial mean
P1 <- matrix(0, 2, 2)       # Initial state covariance (zero for diffuse initialization)
P1inf <- diag(2)            # Diffuse indicator matrix (identity matrix for both states)

# Define the state-space model
model <- SSModel(
  y ~ -1 + SSMcustom(
    Z = Z_array,    # Observation matrix: y(t) = beta(t)*x(t) + noise
    T = T_matrix,   # State transition matrix
    R = diag(2),    # State noise selection matrix
    Q = Q_matrix,   # State noise covariance matrix
    a1 = a1,        # Initial mean
    P1 = P1,        # Initial covariance
    P1inf = P1inf   # Selection matrix for the exact diffuse initialisation 
  ),
  H = NA              # Observation noise variance to be estimated
)

# Estimate the parameters via MLE
fit <- fitSSM(model, inits=c(0.1, 0.1), method="BFGS")

sqrt(fit$model$H[1, 1, 1])
sqrt(fit$model$Q[2, 2, 1])
# Get filtered estimates (the smoother is not needed here, as I only need the last value)
filtered <- KFS(fit$model, smoothing="none")
ts.plot(cbind(beta, filtered$att[,1]),col=1:2)

You might to test with different initial values in case you are getting trapped to some poor local optimum.

[fipelle]

Hi,

I have already tried a few different initial values with no luck. I suppose that my poor optimization results might be due to either the BFGS algorithm (or its tolerance) or the fact that I am estimating both H and the last element of Q. I will try fixing the last element of Q to a small value and estimate H by targeting the signal-to-noise ratio with the optimizer. If I want to change the tolerance, how do I do it? Alternatively, can I constraint H to be bigger than the last element of Q?

[helske]

Well as you can see from my example above, the model should work ok if the data is actually generated according to the model specification. Are you sure you do not need intercept in your model?

fitSSM uses optim in the background, so you can switch to different algorithm and control tolerances etc of optim by passing them as arguments to fitSSM. There are also some examples in the docs of fitSSM how to create own model update functions which can used for example to build the constraints you are thinking. Alternatively, you can build your own objective function by using logLik method and then you can use arbitrary optimization routines available in R. There is an example of this in ?KFAS.