In josephsdavid/tstest: tstest

Forecasting

Some helpers:

tswgen <- function(n,sn=0){
    sig <- function(...){
        gen.sigplusnoise.wge(n=n,...,sn=sn)
    }
    ari <- function(...,dif){
        gen.arima.wge(n=n,...,sn=sn)
    }
    aru <- function(...){
        gen.aruma.wge(n=n,..., sn=sn)
    }

    arma <- function(...){
        gen.arma.wge(n=n,..., sn=sn)
    }
    list("sig"=sig,"ari"=ari,"aru"=aru, "arma" = arma)
}

Forecasts from Signal-Plus-Noise

$$X_t = s_t + Z_t$$ Where:

• $s_t$ is a deterministic signal
• $Z_t$ is a zero-mean, stationary process

Examples:

• $s_t = a + b t$ is a linear signal

• $s_t = a + b t + c t^2$ is a quadratic signal

• $s_t = \cos \left(2 \pi f t + C \right)$

Forecasting setting

Forecast predicts future behavior given finite reaizations. ARMA forecasting has gotten really populat lately.

• Curve fitting

  • Has the underlying assumption that future behavior follows a deterministic path

• Time series based forecasting (ARMA et al)

  • Underlying assumption that the future is only guided by its correlation to the past.

Box-Jenkins approach

We use this

Strategy + Notation

We normally have: $$X_1, X_2, ..., X_{t_0}$$ Our goal is to forecast a future value, say $X_{t_{0+4}}$

Typicall $t_0 = n$ ie it is the last value of the realization. However, sometimes to evaluate our forecast, we compare it to the last few values of the realization.

Notation

$\widehat{X_{t_0}}(\ell)$ is the forecast of $X_{t_{0+\ell}}$ given data up to time $t_0$

• $t_0$ is called the forecast origin

• $\ell$ is called the lead time, the number of steps ahead which we want to forecast

some math

$$X_t - \phi_1 X_{t-1} = (1-\phi_1) \mu + a_t$$

• assume we know $\phi_1$

• we don-t know X11, we dont know mu, and we dont know a12

Iterative forecasting

We assume a in the future is 0, then we have

knitr::read_chunk('R/forecastit.R')

<<forecastfun>>

ex74 <- ar1forcastgen(phi=.8, mu = 24.17)
times <- c(1,2,3,4,5)
lapply(times, ex74, xprev = 22.93) %>% as.data.frame
ex74(22.93,l=2)

Discussion of Ar1 forcasts

AR1 forcasts damp to the sample mean.

AR(p) forecast math

$$\hat{X}{t_0} (\ell) = \phi_1 \hat{X}{t_0}(\ell -1) + ... + \phi_p \hat{X}_{t_0}(\ell - p) + \bar{X} (1 - \phi_1 - ... - \phi_p))$$

AR(p) forecasts:

<<arpfun>>
x <- c(27.7,23.4,21.2,21.1,22.7)
p <- c(1.6,-.8)
avg <- 29.4
arp(phi = p,vec = x,l = 5,mu = avg)

With tswge

positive phi, AR(1)

<<tswarplus>>

Negative phi, AR(1)

<<tswarmin>>

It has the same behavior as the AR(functions). Sweet.

AR(2)

<<tswar2>>

Eventual Forecast Function:

As $\ell$ gets large, it is most reasonable to simply forecast the sample mean, because of how they converge

ARMA(p,q)

Some math

$$\hat{X}{t_0} (\ell) = \sum{j = 1}^p \phi_j \hat{X}{t_0}(\ell - j) - \sum{j = \ell}^q \theta_j \hat{a_{t_0 + \ell - j}} + \bar{X} \left[ 1- \sum_{i = 1}^p \phi_i \right]$$ Where the a hat is incredibly hard to calculate (we have to do back casting, which we can worry about after the test). Challeng: Write this function recursively (is it possible?????) in R. We are just going to use fore.arma.wge for now

ts7524 <- tswgen(75, sn = 24)
p <- c(1.6, -.8)
thet <- (-.9)
x1 <- ts7524$arma(phi = p)
fore.arma.wge(x1, phi = p, n.ahead = 20, limits =F)
x2 <- ts7524$arma(phi = p, theta = thet)
fore.arma.wge(x2, phi = p, theta = thet, n.ahead = 20, limits =F)

Example

data(llynx)
plotts.wge(llynx)

fore.arma.wge(llynx, phi = c(0.7,0.1,-0.2,-0.3),theta = -.6, n.ahead = 20, limits = F)

psi weights

Here is basically a useless page because i am terrible at bookdown

psi.weights.wge(phi = c(0.4,-.6,.8), lag.max = 5)
psiweights(phi = c(0.4,-.6,.8),  l = 6)

Probability limits

First, go back and review the GLP, and understanf that ARMA processes can be expressed in GLP form as:

$$X_t - \mu = \sum_{j = 0}^\infty \psi_j a_{t-j}$$, where at is normal white noise. some more math: the forecast error is defined as:

$$e_{t_0} (\ell) = X_{t_0 + \ell} - \hat{X}{t_0} (\ell)$$ We can rewrite this in GLP form as $$ = \sum{k = 0}^{\ell -1} \psi_k a_{t_0 + \ell - k}$$ Since a is normal, e should also be normal, its expeced value should be zero, and its variance is just the variance of it in GLP form. We can rewrite this as: $$\mathrm{var}(e_{t_0} (\ell)) = \sigma_a^2 \sum_{j = 0}^{\ell -1} \psi_j^2$$ This saves us a bunch of time and we just need to calculate the white noise variance. This is hard so we can just get it from one of the outputs of fore.arma.wge, and then do the math from there.

An alternatve way to calculate psi weights

$$\phi(B)X_t = \theta(B)a_t \rightarrow X_t = \frac{\theta(B)}{\psi(B)}a_t = \sum_{j = 0}^\infty \psi_j B^j a_t$$ Then you just do the long division and get your answers, by dividing off that tricky a_t component.

Calculating them

Since we know all this, we have that $$\frac{e_{t_0}(\ell) - 0}{\left[\sigma_a^2 \sum_{j = 0}^{\ell -1} \psi_j^2\right]^{1/2}} \approx N(0,1)$$ $$=\frac{X_{t_0 + \ell} - \hat{X}{t_0} (\ell)}{ {\left[\sigma_a^2 \sum{j = 0}^{\ell -1} \psi_j^2\right]^{1/2}} }$$

Given the above, we have that the probability of the above lying in betwieen $\pm 1.96$ is about 0.95. Then we expand that out, rewriting it/

$$\implies P\left(-1.96 \leq \frac{X_{t_0 + \ell} - \hat{X}{t_0} (\ell)}{ {\left[\sigma_a^2 \sum{j = 0}^{\ell -1} \psi_j^2\right]^{1/2}} \leq 1.96}\right)$$ $$ \implies P\left(-1.96 \left[\sigma_a^2 \sum_{j = 0}^{\ell -1} \psi_j^2\right]^{1/2} \leq X_{t_0 + \ell} - \hat{X}{t_0} (\ell) \leq 1.96\left[\sigma_a^2 \sum{j = 0}^{\ell -1} \psi_j^2\right]^{1/2} \right) = 0.95$$

We then have

$$ \implies P\left(\hat{X}{t_0} (\ell)-1.96 \left[\sigma_a^2 \sum{j = 0}^{\ell -1} \psi_j^2\right]^{1/2} \leq X_{t_0 + \ell} \leq \hat{X}{t_0} (\ell) +1.96\left[\sigma_a^2 \sum{j = 0}^{\ell -1} \psi_j^2\right]^{1/2} \right) = 0.95$$

This gives us our result!!! Note that this is a huge simplification of the actual method of doing this, it is actually much more complicated what our computer does.

ASE

<<tswase>>
res1 <- tswase(vec = llynx, phi = c(1.3,-0.7,0.1,-0.2), ahead = 12, lastn =T, limits = F)
library(pander)
res2 <- tswase(vec = llynx, phi = c(0.7,0.1,-0.2,-0.3),theta = -.6, ahead = 12, lastn =T, limits = F)


pander(res1)
pander(res2)

ARIMA forecasts

These are too hard to do by hand, we do not need to worry about it. We will use fore.aruma/fore.arima instead. However, it is interesting as an excercise to do an ARIMA(0,1,0) forecast

An example:

Consider $$(1-0.9B)C_t = a_t$$ We have that $\phi_1$ = 1, as is the case with these models. If we recall the general AR(p) forecast equation:

$$\hat{X}{t_0} (\ell) = \phi_1 \hat{X}{t_0}(\ell -1) + ... + \phi_p \hat{X}{t_0}(\ell - p) + \bar{X} (1 - \phi_1 - ... - \phi_p))$$ We see that the right hand term with Xbar goes to zero, and the phi value is one. This gives us $$\hat{X}{t_0} (\ell) = \hat{X}_{t_0} (\ell - 1)$$

Which says that we are just going to forecast the last value we obseved over and over and over. Pretty boring, but that is fine with me. Note that the probability limits of these forecasts increase in width with no bound. Pretty crap forecast huh.

Seasonal Forecasts

Consider a quarterly model:

$$(1 - B^4) (X_t - \mu) = a_t$$

If we rewrite the general AR forecast, since we still have a root at 1 (a fourth root), the xbar component still goes away.Then, we are left with

$$\hat{X}{t_0} (\ell) = \hat{X}{t_0} (\ell - 4)$$

This again makes a lot of sense. Each forecast is exactly the same as it was 4 times ago. A pretty simple and straightforward deal. We will after the test combine these models.

Airline models

In order to forecast into the future with a seasonal model, we want to not only include the existing seasonal trend, but also a 1-B wandering trend. Seasonal and ARIMA components. This way we get forecasts going off of themselves. This works because of

factor.wge(phi = c(rep(0,11),1))

We see that we have a 1-B component in there already, which means we can put it in there and it will get cancelled out without messing up out model!!!

The most general seasonal model we can have is:

$$\phi(B)(1-B)(1-B^s)(X_t - \mu) = \theta(B)a_t$$

The forecasts follow the example of the last s points, but it is not a replication. Models like this are known as airline models, specifically with s = 12. This is how you forecast the trend and forecast the seasonality.

signal plus noise forecasts

Linear Signal

• Fit $b_0$ and $b_1$ using least squares (or true values)

• Find residuals

• Fit an AR(p) to the residuals

Then find forecasts on the next noise value, generating an estimate of the signal (with the regression) and the noise (with the residuals). We will get more into this later and use some psis to find the limits.

Cyclic signal

Will do later, for now for both we can just do fore.sigplusnoise.wge For early lags, the noise behavior is dominant, and as we get further out, the signal dominates

josephsdavid/tstest documentation built on July 3, 2019, 4:48 p.m.