Title: | Improved Prediction Intervals for ARIMA Processes and Structural Time Series |
---|---|
Description: | Prediction intervals for ARIMA and structural time series models using importance sampling approach with uninformative priors for model parameters, leading to more accurate coverage probabilities in frequentist sense. Instead of sampling the future observations and hidden states of the state space representation of the model, only model parameters are sampled, and the method is based solving the equations corresponding to the conditional coverage probability of the prediction intervals. This makes method relatively fast compared to for example MCMC methods, and standard errors of prediction limits can also be computed straightforwardly. |
Authors: | Jouni Helske |
Maintainer: | Jouni Helske <[email protected]> |
License: | GPL-3 |
Version: | 1.0.4 |
Built: | 2024-11-14 03:43:00 UTC |
Source: | https://github.com/helske/tsPI |
Function acv_arma
computes a theoretical autocovariance function of ARMA process.
acv_arma(phi, theta, n)
acv_arma(phi, theta, n)
phi |
vector containing the AR parameters |
theta |
vector containing the MA parameters |
n |
length of the time series |
vector of length n containing the autocovariances
## Example from Brockwell & Davis (1991, page 92-94) ## also in help page of ARMAacf (from stats) n <- 0:9 answer <- 2^(-n) * (32/3 + 8 * n) /(32/3) acv <- acv_arma(c(1.0, -0.25), 1.0, 10) all.equal(acv/acv[1], answer)
## Example from Brockwell & Davis (1991, page 92-94) ## also in help page of ARMAacf (from stats) n <- 0:9 answer <- 2^(-n) * (32/3 + 8 * n) /(32/3) acv <- acv_arma(c(1.0, -0.25), 1.0, 10) all.equal(acv/acv[1], answer)
Function arima_pi
computes prediction intervals for ARIMA processes
with exogenous variables using importance sampling. For regression coefficients,
diffuse (uninformative) prior is used, whereas multiple options for
prior distributions for ARMA coefficients are supported.
arima_pi( x, order, xreg = NULL, n_ahead = 1, level = 0.95, median = TRUE, se_limits = TRUE, prior = "uniform", custom_prior, custom_prior_args = NULL, nsim = 1000, invertibility = FALSE, last_only = FALSE, return_weights = FALSE, ... )
arima_pi( x, order, xreg = NULL, n_ahead = 1, level = 0.95, median = TRUE, se_limits = TRUE, prior = "uniform", custom_prior, custom_prior_args = NULL, nsim = 1000, invertibility = FALSE, last_only = FALSE, return_weights = FALSE, ... )
x |
vector containing the time series |
order |
vector of length 3 with values p,d,q corresponding to the number of AR parameters, degree of differencing and number of MA parameters. |
xreg |
matrix or data frame containing the exogenous variables (not including the intercept which is always included for non-differenced series) |
n_ahead |
length of the forecast horizon. |
level |
desired frequentist coverage probability of the prediction intervals. |
median |
compute the median of the prediction interval. |
se_limits |
compute the standard errors of the prediction interval limits. |
prior |
prior to be used in importance sampling for AR and MA parameters.
Defaults to uniform prior. Several Jeffreys' priors are also available (see |
custom_prior |
function for computing custom prior. First argument must be a vector containing the AR and MA parameters (in that order). |
custom_prior_args |
list containing additional arguments to |
nsim |
number of simulations used in importance sampling. Default is 1000. |
invertibility |
Logical, should the priors include invertibility constraint? Default is |
last_only |
compute the prediction intervals only for the last prediction step. |
return_weights |
Return (scaled) weights used in importance sampling. |
... |
additional arguments for |
a list containing the prediction intervals. @references
Helske, J. and Nyblom, J. (2015). Improved frequentist prediction intervals for autoregressive models by simulation. In Siem Jan Koopman and Neil Shephard, editors, Unobserved Components and Time Series Econometrics. Oxford University Press. https://urn.fi/URN:NBN:fi:jyu-201603141839
Helske, J. and Nyblom, J. (2014). Improved frequentist prediction intervals for ARMA models by simulation. In Johan Knif and Bernd Pape, editors, Contributions to Mathematics, Statistics, Econometrics, and Finance: essays in honour of professor Seppo Pynnönen, number 296 in Acta Wasaensia, pages 71–86. University of Vaasa. https://urn.fi/URN:NBN:fi:jyu-201603141836
set.seed(123) x <- arima.sim(n = 30, model = list(ar = 0.9)) pred_arima <- predict(arima(x, order = c(1,0,0)), n.ahead = 10, se.fit = TRUE) pred_arima <- cbind(pred = pred_arima$pred, lwr = pred_arima$pred - qnorm(0.975)*pred_arima$se, upr = pred_arima$pred + qnorm(0.975)*pred_arima$se) pred <- arima_pi(x, order = c(1,0,0), n_ahead = 10) ts.plot(ts.union(x,pred_arima, pred[,1:3]), col = c(1,2,2,2,3,3,3), lty = c(1,1,2,2,1,2,2))
set.seed(123) x <- arima.sim(n = 30, model = list(ar = 0.9)) pred_arima <- predict(arima(x, order = c(1,0,0)), n.ahead = 10, se.fit = TRUE) pred_arima <- cbind(pred = pred_arima$pred, lwr = pred_arima$pred - qnorm(0.975)*pred_arima$se, upr = pred_arima$pred + qnorm(0.975)*pred_arima$se) pred <- arima_pi(x, order = c(1,0,0), n_ahead = 10) ts.plot(ts.union(x,pred_arima, pred[,1:3]), col = c(1,2,2,2,3,3,3), lty = c(1,1,2,2,1,2,2))
arima_pi
Computes expected coverage probabilities of the prediction intervals of ARMA process by simulating time series from the known model.
avg_coverage_arima( phi = NULL, theta = NULL, d = 0, n, n_ahead = 1, nsim2, nsim = 100, level = 0.95, prior = "uniform", return_all_coverages = FALSE, ... )
avg_coverage_arima( phi = NULL, theta = NULL, d = 0, n, n_ahead = 1, nsim2, nsim = 100, level = 0.95, prior = "uniform", return_all_coverages = FALSE, ... )
phi |
vector containing the AR parameters |
theta |
vector containing the MA parameters |
d |
degree of differencing |
n |
length of the time series |
n_ahead |
length of the forecast horizon |
nsim2 |
number of simulations used in computing the expected coverage |
nsim |
number of simulations used in importance sampling |
level |
desired coverage probability of the prediction intervals |
prior |
prior to be used in importance sampling. Multiple choices are allowed. |
return_all_coverages |
return raw results i.e. coverages for each simulations. When |
... |
additional arguments to |
a list containing the coverage probabilities
## Not run: set.seed(123) # takes a while, notice se, increase nsim2 to get more accurate results avg_coverage_arima(phi = 0.9, n = 50, n_ahead = 10, nsim2 = 100) avg_coverage_arima(phi = 0.9, theta = -0.6, n = 50, n_ahead = 10, nsim2 = 100) ## End(Not run)
## Not run: set.seed(123) # takes a while, notice se, increase nsim2 to get more accurate results avg_coverage_arima(phi = 0.9, n = 50, n_ahead = 10, nsim2 = 100) avg_coverage_arima(phi = 0.9, theta = -0.6, n = 50, n_ahead = 10, nsim2 = 100) ## End(Not run)
struct_pi
and plug-in methodComputes expected coverage probabilities of the prediction intervals of structural time series model. Note that for the plug-in method only standard deviations are assumed to be identical to their estimates, but the initial values for the states are still treated as diffuse. Because of this, plug-in method often performs relatively well in case of structural time series models compared to similar type of ARIMA models (local level and local linear trend models are closely related to ARIMA(0,1,1) and ARIMA(0,2,2) models), and in some cases even outperforms the importance sampling approach with uniform prior (see examples). This is not suprising, as local level and local linear trend models are closely related to ARIMA(0,1,1) and ARIMA(0,2,2) models, and the effect of uncertainty in MA components is not as significant as the uncertainty of AR components
avg_coverage_struct( type = c("level", "trend", "BSM"), sds, frequency = 1, n, n_ahead = 1, nsim2, nsim = 100, level = 0.95, prior = "uniform", return_all_coverages = FALSE, ... )
avg_coverage_struct( type = c("level", "trend", "BSM"), sds, frequency = 1, n, n_ahead = 1, nsim2, nsim = 100, level = 0.95, prior = "uniform", return_all_coverages = FALSE, ... )
type |
Type of model. See |
sds |
vector containing the standard deviations of the model (observation error, level, slope, and seasonal). |
frequency |
frequency of the series, needed for seasonal component. |
n |
length of the time series |
n_ahead |
length of the forecast horizon |
nsim2 |
number of simulations used in computing the expected coverage |
nsim |
number of simulations used in importance sampling |
level |
desired coverage probability of the prediction intervals |
prior |
prior to be used in importance sampling. |
return_all_coverages |
return raw results i.e. coverages for each simulations. When |
... |
additional arguments to |
a list containing the coverage probabilities
## Not run: set.seed(123) # takes a while, notice se, increase nsim2 to get more accurate results avg_coverage_struct(type = "level", sds = c(1, 0.1), n = 50, n_ahead = 10, nsim2 = 100) avg_coverage_struct(type = "BSM", sds = c(1, 1, 0.1, 10), frequency = 4, n = 50, n_ahead = 10, nsim2 = 100) ## End(Not run)
## Not run: set.seed(123) # takes a while, notice se, increase nsim2 to get more accurate results avg_coverage_struct(type = "level", sds = c(1, 0.1), n = 50, n_ahead = 10, nsim2 = 100) avg_coverage_struct(type = "BSM", sds = c(1, 1, 0.1, 10), frequency = 4, n = 50, n_ahead = 10, nsim2 = 100) ## End(Not run)
Function dacv_arma
computes the partial derivatives of theoretical autocovariance function of ARMA process
dacv_arma(phi, theta, n)
dacv_arma(phi, theta, n)
phi |
vector containing the AR parameters |
theta |
vector containing the MA parameters |
n |
length of the time series |
matrix containing the partial derivatives autocovariances, each column corresponding to one parameter of vector (phi,theta) (in that order)
Fortran implementation of InformationMatrixARMA
function of
FitARMA
package, except that the function uses the same
ARMA model definition as arima
, where both the
AR and MA parts of the model are on the right side of the equation, i.e.
MA coefficients differ in sign compared to InformationMatrixARMA
.
information_arma(phi = NULL, theta = NULL)
information_arma(phi = NULL, theta = NULL)
phi |
Autoregressive coefficients. |
theta |
Moving average coefficients. |
Large sample approximation of information matrix for ARMA process.
Box, G. and Jenkins, G. (1970). Time Series Analysis: Forecasting and Control. San Francisco: Holden-Day.
McLeod, A. I. and Zhang, Y., (2007). Faster ARMA maximum likelihood estimation Computational Statistics & Data Analysis 52(4) URL https://dx.doi.org/10.1016/j.csda.2007.07.020
These functions compute different types of importance weights based on Jeffreys's priors used in arima_pi
.
approx_joint_jeffreys(psi, xreg = NULL, p, q, n) approx_marginal_jeffreys(psi, p, q) exact_joint_jeffreys(psi, xreg = NULL, p, q, n) exact_marginal_jeffreys(psi, p, q, n)
approx_joint_jeffreys(psi, xreg = NULL, p, q, n) approx_marginal_jeffreys(psi, p, q) exact_joint_jeffreys(psi, xreg = NULL, p, q, n) exact_marginal_jeffreys(psi, p, q, n)
psi |
vector containing the ar and ma parameters (in that order). |
xreg |
matrix or data frame containing the exogenous variables (not including the intercept which is always included for non-differenced series) |
p |
number of ar parameters |
q |
number of ma parameters |
n |
length of the time series |
Function struct_pi
computes prediction intervals for structural time series
with exogenous variables using importance sampling.
struct_pi( x, type = c("level", "trend", "BSM"), xreg = NULL, n_ahead = 1, level = 0.95, median = TRUE, se_limits = TRUE, prior = "uniform", custom_prior, custom_prior_args = NULL, nsim = 1000, inits = NULL, last_only = FALSE, return_weights = FALSE )
struct_pi( x, type = c("level", "trend", "BSM"), xreg = NULL, n_ahead = 1, level = 0.95, median = TRUE, se_limits = TRUE, prior = "uniform", custom_prior, custom_prior_args = NULL, nsim = 1000, inits = NULL, last_only = FALSE, return_weights = FALSE )
x |
vector containing the time series |
type |
type of model. Possible options are |
xreg |
matrix or data frame containing the exogenous variables (not including the intercept which is always included for non-differenced series) |
n_ahead |
length of the forecast horizon. |
level |
desired frequentist coverage probability of the prediction intervals. |
median |
compute the median of the prediction interval. |
se_limits |
compute the standard errors of the prediction interval limits. |
prior |
prior to be used in importance sampling for log-sd parameters. Defaults to uniform prior on logarithm of standard deviations (with constraints that all variances are smaller than 1e7). If "custom", a user-defined custom prior is used (see next arguments). |
custom_prior |
function for computing custom prior. First argument must be a vector containing the log-variance parameters (observation error, level, slope, and seasonal). |
custom_prior_args |
list containing additional arguments to |
nsim |
number of simulations used in importance sampling. Default is 1000. |
inits |
initial values for log-sds |
last_only |
compute the prediction intervals only for the last prediction step. |
return_weights |
Return (scaled) weights used in importance sampling. |
a list containing the prediction intervals.
Helske, J. (2015). Prediction and interpolation of time series by state space models. University of Jyväskylä. PhD thesis, Report 152. https://urn.fi/URN:NBN:fi:jyu-201603111829
pred_StructTS <- predict(StructTS(Nile, type ="level"), n.ahead = 10, se.fit = TRUE) pred_StructTS <- cbind(pred = pred_StructTS$pred, lwr = pred_StructTS$pred - qnorm(0.975)*pred_StructTS$se, upr = pred_StructTS$pred + qnorm(0.975)*pred_StructTS$se) set.seed(123) pred <- struct_pi(Nile, type = "level", n_ahead = 10) ts.plot(ts.union(Nile,pred_StructTS, pred[,1:3]), col = c(1,2,2,2,3,3,3), lty = c(1,1,2,2,1,2,2))
pred_StructTS <- predict(StructTS(Nile, type ="level"), n.ahead = 10, se.fit = TRUE) pred_StructTS <- cbind(pred = pred_StructTS$pred, lwr = pred_StructTS$pred - qnorm(0.975)*pred_StructTS$se, upr = pred_StructTS$pred + qnorm(0.975)*pred_StructTS$se) set.seed(123) pred <- struct_pi(Nile, type = "level", n_ahead = 10) ts.plot(ts.union(Nile,pred_StructTS, pred[,1:3]), col = c(1,2,2,2,3,3,3), lty = c(1,1,2,2,1,2,2))
Package tsPI computes prediction intervals for ARIMA and structural time series models by using importance sampling approach with uninformative priors for model parameters, leading to more accurate coverage probabilities in frequentist sense. Instead of sampling the future observations and hidden states of the state space representation of the model, only model parameters are sampled, and the method is based solving the equations corresponding to the conditional coverage probability of the prediction intervals. This makes method relatively fast compared to for example MCMC methods, and standard errors of prediction limits can also be computed straightforwardly.
Helske, J. and Nyblom, J. (2013). Improved frequentist prediction intervals for autoregressive models by simulation. In Siem Jan Koopman and Neil Shephard, editors, Unobserved Components and Time Series Econometrics. Oxford University Press. In press.
Helske, J. and Nyblom, J. (2014). Improved frequentist prediction intervals for ARMA models by simulation. In Johan Knif and Bernd Pape, editors, Contributions to Mathematics, Statistics, Econometrics, and Finance: essays in honour of professor Seppo Pynnönen, number 296 in Acta Wasaensia, pages 71–86. University of Vaasa.
Helske, J. (2015). Prediction and interpolation of time series by state space models. University of Jyväskylä. PhD thesis, Report 152.