Description Usage Arguments Details Value Examples

Construct an object of class `ssm_ulg`

by directly defining the corresponding terms of
the model.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 | ```
ssm_ulg(
y,
Z,
H,
T,
R,
a1,
P1,
init_theta = numeric(0),
D,
C,
state_names,
update_fn = default_update_fn,
prior_fn = default_prior_fn
)
``` |

`y` |
Observations as time series (or vector) of length |

`Z` |
System matrix Z of the observation equation as m x 1 or m x n matrix. |

`H` |
Vector of standard deviations. Either a scalar or a vector of length n. |

`T` |
System matrix T of the state equation. Either a m x m matrix or a m x m x n array. |

`R` |
Lower triangular matrix R the state equation. Either a m x k matrix or a m x k x n array. |

`a1` |
Prior mean for the initial state as a vector of length m. |

`P1` |
Prior covariance matrix for the initial state as m x m matrix. |

`init_theta` |
Initial values for the unknown hyperparameters theta. |

`D` |
Intercept terms for observation equation, given as a length n vector. |

`C` |
Intercept terms for state equation, given as m x n matrix. |

`state_names` |
Names for the states. |

`update_fn` |
Function which returns list of updated model components given input vector theta. See details. |

`prior_fn` |
Function which returns log of prior density given input vector theta. |

The general univariate linear-Gaussian model is defined using the following observational and state equations:

*y_t = D_t + Z_t α_t + H_t ε_t, (\textrm{observation equation})*

*α_{t+1} = C_t + T_t α_t + R_t η_t, (\textrm{transition equation})*

where *ε_t \sim N(0, 1)*, *η_t \sim N(0, I_k)* and
*α_1 \sim N(a_1, P_1)* independently of each other.
Here k is the number of disturbance terms which can be less than m, the number of states.

The `update_fn`

function should take only one
vector argument which is used to create list with elements named as
`Z`

, `H`

`T`

, `R`

, `a1`

, `P1`

, `D`

, and `C`

,
where each element matches the dimensions of the original model.
If any of these components is missing, it is assumed to be constant wrt. theta.
Note that while you can input say R as m x k matrix for `ssm_ulg`

,
`update_fn`

should return R as m x k x 1 in this case.
It might be useful to first construct the model without updating function and then check
the expected structure of the model components from the output.

Object of class `ssm_ulg`

.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 | ```
# Regression model with time-varying coefficients
set.seed(1)
n <- 100
x1 <- rnorm(n)
x2 <- rnorm(n)
b1 <- 1 + cumsum(rnorm(n, sd = 0.5))
b2 <- 2 + cumsum(rnorm(n, sd = 0.1))
y <- 1 + b1 * x1 + b2 * x2 + rnorm(n, sd = 0.1)
Z <- rbind(1, x1, x2)
H <- 0.1
T <- diag(3)
R <- diag(c(0, 1, 0.1))
a1 <- rep(0, 3)
P1 <- diag(10, 3)
# updates the model given the current values of the parameters
update_fn <- function(theta) {
R <- diag(c(0, theta[1], theta[2]))
dim(R) <- c(3, 3, 1)
list(R = R, H = theta[3])
}
# prior for standard deviations as half-normal(1)
prior_fn <- function(theta) {
if(any(theta < 0)){
log_p <- -Inf
} else {
log_p <- sum(dnorm(theta, 0, 1, log = TRUE))
}
log_p
}
model <- ssm_ulg(y, Z, H, T, R, a1, P1,
init_theta = c(1, 0.1, 0.1),
update_fn = update_fn, prior_fn = prior_fn)
out <- run_mcmc(model, iter = 10000)
out
sumr <- summary(out, variable = "state")
ts.plot(sumr$Mean, col = 1:3)
lines(b1, col= 2, lty = 2)
lines(b2, col= 3, lty = 2)
# Perhaps easiest way to construct a general SSM for bssm is to use the
# model building functionality of KFAS:
library("KFAS")
model_kfas <- SSModel(log(drivers) ~ SSMtrend(1, Q = 5e-4)+
SSMseasonal(period = 12, sea.type = "trigonometric", Q = 0) +
log(PetrolPrice) + law, data = Seatbelts, H = 0.005)
# use as_bssm function for conversion, kappa defines the
# prior variance for diffuse states
model_bssm <- as_bssm(model_kfas, kappa = 100)
# define updating function for parameter estimation
# we can use SSModel and as_bssm functions here as well
# (for large model it is more efficient to do this
# "manually" by constructing only necessary matrices,
# i.e., in this case a list with H and Q)
updatefn <- function(theta){
model_kfas <- SSModel(log(drivers) ~ SSMtrend(1, Q = theta[1]^2)+
SSMseasonal(period = 12,
sea.type = "trigonometric", Q = theta[2]^2) +
log(PetrolPrice) + law, data = Seatbelts, H = theta[3]^2)
as_bssm(model_kfas, kappa = 100)
}
prior <- function(theta) {
if(any(theta < 0)) -Inf else sum(dnorm(theta, 0, 0.1, log = TRUE))
}
init_theta <- rep(1e-2, 3)
c("sd_level", "sd_seasonal", "sd_y")
model_bssm <- as_bssm(model_kfas, kappa = 100,
init_theta = init_theta,
prior_fn = prior, update_fn = updatefn)
## Not run:
out <- run_mcmc(model_bssm, iter = 10000, burnin = 5000)
out
# Above the regression coefficients are modelled as time-invariant latent states.
# Here is an alternative way where we use variable D so that the
# coefficients are part of parameter vector theta:
updatefn2 <- function(theta) {
# note no PetrolPrice or law variables here
model_kfas2 <- SSModel(log(drivers) ~ SSMtrend(1, Q = theta[1]^2)+
SSMseasonal(period = 12, sea.type = "trigonometric", Q = theta[2]^2),
data = Seatbelts, H = theta[3]^2)
X <- model.matrix(~ -1 + law + log(PetrolPrice), data = Seatbelts)
D <- t(X %*% theta[4:5])
as_bssm(model_kfas2, D = D, kappa = 100)
}
prior2 <- function(theta) {
if(any(theta[1:3] < 0)) {
-Inf
} else {
sum(dnorm(theta[1:3], 0, 0.1, log = TRUE)) +
sum(dnorm(theta[4:5], 0, 10, log = TRUE))
}
}
init_theta <- c(rep(1e-2, 3), 0, 0)
names(init_theta) <- c("sd_level", "sd_seasonal", "sd_y", "law", "Petrol")
model_bssm2 <- updatefn2(init_theta)
model_bssm2$theta <- init_theta
model_bssm2$prior_fn <- prior2
model_bssm2$update_fn <- updatefn2
out2 <- run_mcmc(model_bssm2, iter = 10000, burnin = 5000)
out2
## End(Not run)
``` |

Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.