Forecasting of Stationary and Non-Stationary Time Series

Description

Methods to compute linear h-step ahead prediction coefficients based on localised and iterated Yule-Walker estimates and empirical mean squared and absolute prediction errors for the resulting predictors. Also, functions to compute autocovariances for AR(p) processes, to simulate tvARMA(p,q) time series, and to verify an assumption from Kley et al. (2019).

Details

Contents

The core functionality of this R package is accessable via the function predCoef, which is used to compute the linear prediction coefficients, and the functions MSPE and MAPE, which are used to compute the empirical mean squared or absolute prediction errors. Further, the function f can be used to verify condition (10) of Theorem 3.1 in Kley et al. (2019) for any given tvAR(p) model. The function tvARMA can be used to simulate time-varying ARMA(p,q) time series. The function acfARp computes the autocovariances of a AR(p) process from the coefficients and innovations standard deviation.

Author(s)

Tobias Kley

References

Kley, T., Preuss, P. & Fryzlewicz, P. (2019). Predictive, finite-sample model choice for time series under stationarity and non-stationarity. Electronic Journal of Statistics, forthcoming. [cf. https://arxiv.org/abs/1611.04460]

Compute autocovariances of an AR(p) process

Description

This functions returns the autocovariances Cov(X_{t-k}, X_t) of a stationary time series (Y_t) that fulfills the following equation:

Y_t = \sum_{j=1}^p a_j Y_{t-j} + \sigma \varepsilon_{t},

where \sigma > 0, \varepsilon_t is white noise and a_1, \ldots, a_p are real numbers satisfying that the roots z_0 of the polynomial 1 - \sum_{j=1}^p a_j z^j lie strictly outside the unit circle.

Usage

acfARp(a = NULL, sigma, k)

Arguments

Value

Returns autocovariance at lag k of the AR(p) process.

Examples

## Taken from Section 6 in Dahlhaus (1997, AoS)
a1 <- function(u) {1.8 * cos(1.5 - cos(4*pi*u))}
a2 <- function(u) {-0.81}
# local autocovariance for u === 1/2: lag 1
acfARp(a = c(a1(1/2), a2(1/2)), sigma = 1, k = 1)
# local autocovariance for u === 1/2: lag -2
acfARp(a = c(a1(1/2), a2(1/2)), sigma = 1, k = -1)
# local autocovariance for u === 1/2: the variance
acfARp(a = c(a1(1/2), a2(1/2)), sigma = 1, k = 0)

Mean Squared Prediction Errors, for a single h

Description

This function computes the estimated mean squared prediction errors from a given time series and prediction coefficients

Arguments

Details

The array of prediction coefficients coef is expected to be of dimension P x P x H x length(N) x length(t) and in the format as it is returned by the function predCoef. More precisely, for p=1,\ldots,P and the j.Nth element of N element of N the coefficient of the h-step ahead predictor for X_{i+h} which is computed from the observations X_i, \ldots, X_{i-p+1} has to be available via coef[p, 1:p, h, j.N, t==i].

Note that t have to be the indices corresponding to the coefficients.

The resulting mean squared prediction error

\frac{1}{|\mathcal{T}|} \sum_{t \in \mathcal{T}} (X_{t+h} - (X_t, \ldots, X_{t-p+1}) \hat v_{N[j.N],T}^{(p,h)}(t))^2

is then stored in the resulting matrix at position (p, j.N).

Value

Returns a P x length(N) matrix with the results.

Compute f(\delta) for a tvAR(p) process

Description

This functions computes the quantity f(\delta) defined in (24) of Kley et al. (2019) when the underlying process follows an tvAR(p) process. Recall that, to apply Theorem 3.1 in Kley et al. (2019), the function f(\delta) is required to be positive, which can be verified with the numbers returned from this function. The function returns a vector with elements f(\delta) for each \delta in which.deltas, with f(\delta) defined as

f(\delta) := \min_{p_1,p_2 = 0, \ldots, p_{\max}} \min_{N \in \mathcal{N}} \Big| {\rm MSPE}_{s_1/T,m/T}^{(p_1,h)}(\frac{s_1}{T}) - (1+\delta) \cdot {\rm MSPE}_{N/T,m/T}^{(p_2,h)}(\frac{s_1}{T}) \Big|, \quad \delta \geq 0

where T, m, p_{\max}, h are positive integers, \mathcal{N} \subset \{p_{\max}+1, \ldots, T-m-h\}, and s_1 := T-m-h+1.

Usage

f(which.deltas, p_max, h, T, Ns, m, a, sigma)

Arguments

Details

The function {\rm MSPE}_{\Delta_1, \Delta_2}^{(p,h)}(u) is defined, for real-valued u and \Delta_1, \Delta_2 \geq 0, in terms of the second order properties of the process:

{\rm MSPE}_{\Delta_1, \Delta_2}^{(p,h)}(u) := \int_0^1 g^{(p,h)}_{\Delta_1}\Big( u + \Delta_2 (1-x) \Big) {\rm d}x,

with g^{(0,h)}_{\Delta}(u) := \gamma_0(u) and, for p = 1, 2, \ldots,

g^{(p,h)}_{\Delta}(u) := \gamma_0(u) - 2 \big( v_{\Delta}^{(p,h)}(u) \big)' \gamma_0^{(p,h)}(u) + \big( v_{\Delta}^{(p,h)}(u) \big)' \Gamma_0^{(p)}(u) v_{\Delta}^{(p,h)}(u)

\gamma_0^{(p,h)}(u) := \big( \gamma_h(u), \ldots, \gamma_{h+p-1}(u) \big)',

where

v^{(p,h)}_{\Delta}(u) := e'_1 \big( e_1 \big( a_{\Delta}^{(p)}(t) \big)' + H \big)^h,

with e_1 and H defined in the documentation of predCoef and, for every real-valued u and \Delta \geq 0,

a^{(p)}_{\Delta}(u) := \Gamma^{(p)}_{\Delta}(u)^{-1} \gamma^{(p)}_{\Delta}(u),

where

\gamma^{(p)}_{\Delta}(u) := \int_0^1 \gamma^{(p)}(u+\Delta (x-1)) {\rm d}x, \quad \gamma^{(p)}(u) := [\gamma_1(u)\;\ldots\;\gamma_p(u)]',

\Gamma^{(p)}_{\Delta}(u) := \int_0^1 \Gamma^{(p)}(u+\Delta (x-1)) {\rm d}x, \quad \Gamma^{(p)}(u) := (\gamma_{i-j}(u);\,i,j=1,\ldots,p).

The local autocovariances \gamma_k(u) are defined as the lag-k autocovariances of an AR(p) process which has coefficients a_1(u), \ldots, a_p(u) and innovations with variance \sigma(u)^2, because the underlying model is assumed to be tvAR(p)

Y_{t,T} = \sum_{j=1}^p a_j(t/T) Y_{t-j,T} + \sigma(t/T) \varepsilon_{t},

where a_1, \ldots, a_p are real valued functions (defined on [0,1]) and \sigma is a positive function (defined on [0,1]).

Value

Returns a vector with the values f(\delta), as defined in (24) of Kley et al. (2019), where it is now denoted by q(\delta), for each \delta in which.delta.

Examples

## Not run: 
## because computation is quite time-consuming.
n <- 100
a <- list( function(u) {return(0.8+0.19*sin(4*pi*u))} )
sigma <- function (u) {return(1)}

Ns <- seq( floor((n/2)^(4/5)), floor(n^(4/5)),
           ceiling((floor(n^(4/5)) - floor((n/2)^(4/5)))/25) )
which.deltas <- c(0, 0.01, 0.05, 0.1, 0.15, 0.2, 0.4, 0.6)
P_max <- 7
H <- 1
m <- floor(n^(.85)/4)

# now replicate some results from Table 4 in Kley et al. (2019)
f( which.deltas, P_max, h = 1, n - m, Ns, m, a, sigma )
f( which.deltas, P_max, h = 5, n - m, Ns, m, a, sigma )

## End(Not run)

Mean squared or absolute h-step ahead prediction errors

Description

The function MSPE computes the empirical mean squared prediction errors for a collection of h-step ahead, linear predictors (h=1,\ldots,H) of observations X_{t+h}, where m_1 \leq t+h \leq m_2, for two indices m_1 and m_2. The resulting array provides

\frac{1}{m_{\rm lo} - m_{\rm up} + 1} \sum_{t=m_{\rm lo}}^{m_{\rm up}} R_{(t)}^2,

with R_{(t)} being the prediction errors

R_t := | X_{t+h} - (X_t, \ldots, X_{t-p+1}) \hat v_{N,T}^{(p,h)}(t) |,

ordered by magnitude; i.e., they are such that R_{(t)} \leq R_{(t+1)}. The lower and upper limits of the indices are m_{\rm lo} := m_1-h + \lfloor (m_2-m_1+1) \alpha_1 \rfloor and m_{\rm up} := m_2-h - \lfloor (m_2-m_1+1) \alpha_2 \rfloor. The function MAPE computes the empirical mean absolute prediction errors

\frac{1}{m_{\rm lo} - m_{\rm up} + 1} \sum_{t=m_{\rm lo}}^{m_{\rm up}} R_{(t)},

with m_{\rm lo}, m_{\rm up} and R_{(t)} defined as before.

Usage

MSPE(X, predcoef, m1 = length(X)/10, m2 = length(X), P = 1, H = 1,
  N = c(0, seq(P + 1, m1 - H + 1)), trimLo = 0, trimUp = 0)

MAPE(X, predcoef, m1 = length(X)/10, m2 = length(X), P = 1, H = 1,
  N = c(0, seq(P + 1, m1 - H + 1)), trimLo = 0, trimUp = 0)

Arguments

Value

MSPE returns an object of type MSPE that has mspe, an array of size H\timesP\timeslength(N), as an attribute, as well as the parameters N, m1, m2, P, and H. MAPE analogously returns an object of type MAPE that has mape and the same parameters as attributes.

Examples

T <- 1000
X <- rnorm(T)
P <- 5
H <- 1
m <- 20
Nmin <- 20
pcoef <- predCoef(X, P, H, (T - m - H + 1):T, c(0, seq(Nmin, T - m - H, 1)))

mspe <- MSPE(X, pcoef, 991, 1000, 3, 1, c(0, Nmin:(T-m-H)))

plot(mspe, vr = 1, Nmin = Nmin)

Plot a MSPE or MAPE object

Description

The function plot.MSPE plots a MSPE object that is returned by the MSPE function. The function plot.MAPE plots a MAPE object that is returned by the MAPE function.

Usage

## S3 method for class 'MSPE'
plot(x, vr = NULL, h = 1, N_min = 1, legend = TRUE,
  display.mins = TRUE, add.for.legend = 0, ...)

## S3 method for class 'MAPE'
plot(x, vr = NULL, h = 1, N_min = 1, legend = TRUE,
  display.mins = TRUE, add.for.legend = 0, ...)

Arguments

Value

Returns the plot, as specified.

See Also

MSPE, MAPE

h-step Prediction coefficients

Description

This function computes the localised and iterated Yule-Walker coefficients for h-step ahead forecasting of X_{t+h} from X_{t}, ..., X_{t-p+1}, where h = 1, \ldots, H and p = 1, \ldots, P.

Arguments

Details

For every t \in t and every N \in N the (iterated) Yule-Walker estimates \hat v_{N,T}^{(p,h)}(t) are computed. They are defined as

\hat v_{N,T}^{(p,h)}(t) := e'_1 \big( e_1 \big( \hat a_{N,T}^{(p)}(t) \big)' + H \big)^h, \quad N \geq 1,

and

\hat v_{0,T}^{(p,h)}(t) := \hat v_{t,T}^{(p,h)}(t),

with

e_1 := \left(\begin{array}{c} 1 \\ 0 \\ \vdots \\ 0 \end{array} \right), \quad H := \left( \begin{array}{ccccc} 0 & 0 & \cdots & 0 & 0 \\ 1 & 0 & \cdots & 0 & 0 \\ 0 & 1 & \cdots & 0 & 0 \\ \vdots & \ddots & \cdots & 0 & 0 \\ 0 & 0 & \cdots & 1 & 0 \end{array} \right)

and

\hat a_{N,T}^{(p)}(t) := \big( \hat\Gamma_{N,T}^{(p)}(t) \big)^{-1} \hat\gamma_{N,T}^{(p)}(t),

where

\hat\Gamma_{N,T}^{(p)}(t) := \big[ \hat \gamma_{i-j;N,T}(t) \big]_{i,j = 1, \ldots, p}, \quad \hat \gamma_{N,T}^{(p)}(t) := \big( \hat \gamma_{1;N,T}(t), \ldots, \hat \gamma_{p;N,T}(t) \big)'

and

\hat \gamma_{k;N,T}(t) := \frac{1}{N} \sum_{\ell=t-N+|k|+1}^{t} X_{\ell-|k|,T} X_{\ell,T}

is the usual lag-k autocovariance estimator (without mean adjustment), computed from the observations X_{t-N+1}, \ldots, X_{t}.

The Durbin-Levinson Algorithm is used to successively compute the solutions to the Yule-Walker equations (cf. Brockwell/Davis (1991), Proposition 5.2.1). To compute the h-step ahead coefficients we use the recursive relationship

\hat v_{i,N,T}^{(p)}(t,h) = \hat a_{i,N,T}^{(p)}(t) \hat v_{1,N,T}^{(p,h-1)}(t) + \hat v_{i+1,N,T}^{(p,h-1)}(t) I\{i \leq p-1\},

(cf. Section 3.2, Step 3, in Kley et al. (2019)).

Value

Returns a named list with elements coef, t, and N, where coef is an array of dimension P \times P \times H \times length(t) \times length(N), and t, and N are the parameters provided on the call of the function. See the example on how to access the vector \hat v_{N,T}^{(p,h)}(t).

References

Brockwell, P. J. & Davis, R. A. (1991). Time Series: Theory and Methods. Springer, New York.

Examples

T <- 100
X <- rnorm(T)

P <- 5
H <- 1
m <- 20

Nmin <- 25
pcoef <- predCoef(X, P, H, (T - m - H + 1):T, c(0, seq(Nmin, T - m - H, 1)))

## Access the prediction vector for p = 2, h = 1, t = 95, N = 25
p <- 2
h <- 1
t <- 95
N <- 35
res <- pcoef$coef[p, 1:p, h, pcoef$t == t, pcoef$N == N]

Simulation of an tvARMA(p,q) time series.

Description

Returns a simulated time series Y_{1,T}, ..., Y_{T,T} that fulfills the following equation:

Y_{t,T} = \sum_{j=1}^p a_j(t/T) Y_{t-j,T} + \sigma(t/T) \varepsilon_{t} + \sum_{k=1}^q \sigma((t-k)/T) b_k(t/T) \varepsilon_{t-k},

where a_1, \ldots, a_p, b_0, b_1, \ldots, b_q are real-valued functions on [0,1], \sigma is a positive function on [0,1] and \varepsilon_t is white noise.

Usage

tvARMA(T = 128, a = list(), b = list(), sigma = function(u) {    
  return(1) }, innov = function(n) {     rnorm(n, 0, 1) })

Arguments

Value

Returns a tvARMA(p,q) time series with specified parameters.

Examples

## Taken from Section 6 in Dahlhaus (1997, AoS)
a1 <- function(u) {1.8 * cos(1.5 - cos(4 * pi * u))}
a2 <- function(u) {-0.81}
plot(tvARMA(128, a = list(a1, a2), b = list()), type = "l")

Workhorse function for tvARMA time series generation

Description

More explanation!

Arguments

Value

Returns a ...