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) }
$$X_t = s_t + Z_t$$ Where:
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)$
Forecast predicts future behavior given finite reaizations. ARMA forecasting has gotten really populat lately.
Curve fitting
Time series based forecasting (ARMA et al)
We use this
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.
$\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
$$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)
AR1 forcasts damp to the sample mean.
$$\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))$$
<<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)
<<tswarplus>>
<<tswarmin>>
It has the same behavior as the AR(functions). Sweet.
<<tswar2>>
As $\ell$ gets large, it is most reasonable to simply forecast the sample mean, because of how they converge
$$\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)
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)
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)
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.
$$\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.
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.
<<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)
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
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.
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.
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.
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.
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
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.