-
Notifications
You must be signed in to change notification settings - Fork 14
/
kfilter.R
191 lines (177 loc) · 6.94 KB
/
kfilter.R
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
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
#' Kalman Filtering
#'
#' Function \code{kfilter} runs the Kalman filter for the given model,
#' and returns the filtered estimates and one-step-ahead predictions of the
#' states \eqn{\alpha_t} given the data up to time \eqn{t}.
#'
#' For non-Gaussian models, the filtering is based on the approximate
#' Gaussian model.
#'
#' @param model Model of class \code{lineargaussian}, \code{nongaussian} or
#' \code{ssm_nlg}.
#' @param ... Ignored.
#' @return List containing the log-likelihood
#' (approximate in non-Gaussian case), one-step-ahead predictions \code{at}
#' and filtered estimates \code{att} of states, and the corresponding
#' variances \code{Pt} and \code{Ptt} up to the time point n+1 where n is the
#' length of the input time series.
#' @seealso \code{\link{bootstrap_filter}}
#' @export
#' @rdname kfilter
kfilter <- function(model, ...) {
UseMethod("kfilter", model)
}
#' @method kfilter lineargaussian
#' @rdname kfilter
#' @export
#' @examples
#' x <- cumsum(rnorm(20))
#' y <- x + rnorm(20, sd = 0.1)
#' model <- bsm_lg(y, sd_level = 1, sd_y = 0.1)
#' ts.plot(cbind(y, x, kfilter(model)$att), col = 1:3)
kfilter.lineargaussian <- function(model, ...) {
check_missingness(model)
out <- gaussian_kfilter(model, model_type = model_type(model))
colnames(out$at) <- colnames(out$att) <- colnames(out$Pt) <-
colnames(out$Ptt) <- rownames(out$Pt) <- rownames(out$Ptt) <-
names(model$a1)
out$at <- ts(out$at, start = start(model$y), frequency = frequency(model$y))
out$att <- ts(out$att, start = start(model$y), frequency = frequency(model$y))
out
}
#' @method kfilter nongaussian
#' @rdname kfilter
#' @export
kfilter.nongaussian <- function(model, ...) {
kfilter(gaussian_approx(model))
}
#' (Iterated) Extended Kalman Filtering
#'
#' Function \code{ekf} runs the (iterated) extended Kalman filter for the given
#' non-linear Gaussian model of class \code{ssm_nlg},
#' and returns the filtered estimates and one-step-ahead predictions of the
#' states \eqn{\alpha_t} given the data up to time \eqn{t}.
#'
#' @param model Model of class \code{ssm_nlg}.
#' @param iekf_iter Non-negative integer. The default zero corresponds to
#' normal EKF, whereas \code{iekf_iter > 0} corresponds to iterated EKF
#' with \code{iekf_iter} iterations.
#' @return List containing the log-likelihood,
#' one-step-ahead predictions \code{at} and filtered
#' estimates \code{att} of states, and the corresponding variances \code{Pt} and
#' \code{Ptt}.
#' @export
#' @rdname ekf
#' @export
#' @examples
#' \donttest{ # Takes a while on CRAN
#' Sys.setenv("OMP_THREAD_LIMIT" = 2) # For CRAN
#' set.seed(1)
#' mu <- -0.2
#' rho <- 0.7
#' sigma_y <- 0.1
#' sigma_x <- 1
#' x <- numeric(50)
#' x[1] <- rnorm(1, mu, sigma_x / sqrt(1 - rho^2))
#' for(i in 2:length(x)) {
#' x[i] <- rnorm(1, mu * (1 - rho) + rho * x[i - 1], sigma_x)
#' }
#' y <- rnorm(50, exp(x), sigma_y)
#'
#' pntrs <- cpp_example_model("nlg_ar_exp")
#'
#' model_nlg <- ssm_nlg(y = y, a1 = pntrs$a1, P1 = pntrs$P1,
#' Z = pntrs$Z_fn, H = pntrs$H_fn, T = pntrs$T_fn, R = pntrs$R_fn,
#' Z_gn = pntrs$Z_gn, T_gn = pntrs$T_gn,
#' theta = c(mu= mu, rho = rho,
#' log_sigma_x = log(sigma_x), log_sigma_y = log(sigma_y)),
#' log_prior_pdf = pntrs$log_prior_pdf,
#' n_states = 1, n_etas = 1, state_names = "state")
#'
#' out_ekf <- ekf(model_nlg, iekf_iter = 0)
#' out_iekf <- ekf(model_nlg, iekf_iter = 5)
#' ts.plot(cbind(x, out_ekf$att, out_iekf$att), col = 1:3)
#' }
ekf <- function(model, iekf_iter = 0) {
check_missingness(model)
iekf_iter <- check_intmax(iekf_iter, "iekf_iter", positive = FALSE)
out <- ekf_nlg(t(model$y), model$Z, model$H, model$T,
model$R, model$Z_gn, model$T_gn, model$a1, model$P1,
model$theta, model$log_prior_pdf, model$known_params,
model$known_tv_params, model$n_states, model$n_etas,
as.integer(model$time_varying), iekf_iter)
colnames(out$at) <- colnames(out$att) <- colnames(out$Pt) <-
colnames(out$Ptt) <- rownames(out$Pt) <- rownames(out$Ptt) <-
model$state_names
out$at <- ts(out$at, start = start(model$y), frequency = frequency(model$y))
out$att <- ts(out$att, start = start(model$y), frequency = frequency(model$y))
out
}
#' Unscented Kalman Filtering
#'
#' Function \code{ukf} runs the unscented Kalman filter for the given
#' non-linear Gaussian model of class \code{ssm_nlg},
#' and returns the filtered estimates and one-step-ahead predictions of the
#' states \eqn{\alpha_t} given the data up to time \eqn{t}.
#'
#' @param model Model of class \code{ssm_nlg}.
#' @param alpha Positive tuning parameter of the UKF. Default is 0.001. Smaller
#' the value, closer the sigma point are to the mean of the state.
#' @param beta Non-negative tuning parameter of the UKF. The default value is
#' 2, which is optimal for Gaussian states.
#' @param kappa Non-negative tuning parameter of the UKF, which also affects
#' the spread of sigma points. Default value is 0.
#' @return List containing the log-likelihood,
#' one-step-ahead predictions \code{at} and filtered
#' estimates \code{att} of states, and the corresponding variances \code{Pt} and
#' \code{Ptt}.
#' @export
#' @rdname ukf
#' @export
#' @examples
#' \donttest{ # Takes a while on CRAN
#' Sys.setenv("OMP_THREAD_LIMIT" = 2) # For CRAN
#' set.seed(1)
#' mu <- -0.2
#' rho <- 0.7
#' sigma_y <- 0.1
#' sigma_x <- 1
#' x <- numeric(50)
#' x[1] <- rnorm(1, mu, sigma_x / sqrt(1 - rho^2))
#' for(i in 2:length(x)) {
#' x[i] <- rnorm(1, mu * (1 - rho) + rho * x[i - 1], sigma_x)
#' }
#' y <- rnorm(50, exp(x), sigma_y)
#'
#' pntrs <- cpp_example_model("nlg_ar_exp")
#'
#' model_nlg <- ssm_nlg(y = y, a1 = pntrs$a1, P1 = pntrs$P1,
#' Z = pntrs$Z_fn, H = pntrs$H_fn, T = pntrs$T_fn, R = pntrs$R_fn,
#' Z_gn = pntrs$Z_gn, T_gn = pntrs$T_gn,
#' theta = c(mu= mu, rho = rho,
#' log_sigma_x = log(sigma_x), log_sigma_y = log(sigma_y)),
#' log_prior_pdf = pntrs$log_prior_pdf,
#' n_states = 1, n_etas = 1, state_names = "state")
#'
#' out_iekf <- ekf(model_nlg, iekf_iter = 5)
#' out_ukf <- ukf(model_nlg, alpha = 0.01, beta = 2, kappa = 1)
#' ts.plot(cbind(x, out_iekf$att, out_ukf$att), col = 1:3)
#' }
ukf <- function(model, alpha = 0.001, beta = 2, kappa = 0) {
check_missingness(model)
if (alpha <= 0) stop("Parameter 'alpha' should be positive. ")
if (beta < 0) stop("Parameter 'beta' should be non-negative. ")
if (kappa < 0) stop("Parameter 'kappa' should be non-negative. ")
out <- ukf_nlg(t(model$y), model$Z, model$H, model$T,
model$R, model$Z_gn, model$T_gn, model$a1, model$P1,
model$theta, model$log_prior_pdf, model$known_params,
model$known_tv_params, model$n_states, model$n_etas,
as.integer(model$time_varying),
alpha, beta, kappa)
colnames(out$at) <- colnames(out$att) <- colnames(out$Pt) <-
colnames(out$Ptt) <- rownames(out$Pt) <- rownames(out$Ptt) <-
model$state_names
out$at <- ts(out$at, start = start(model$y), frequency = frequency(model$y))
out$att <- ts(out$att, start = start(model$y), frequency = frequency(model$y))
out
}