3
$\begingroup$

I am doing predictions on monthly temperature data for 100 years, from 1901 to 2000 (i.e 1200 data points). I want to know if the method I follow is correct because in my output, I do not see the requisite "randomness" of temperature being reproduced in the prediction.

Here is a link to the plot of the prediction (in red)
https://docs.google.com/file/d/0B1Lm03a_91xiYks5TVJDYU05VUE/edit?usp=sharing

EDIT: added the ACF and PACF of the detrended and de-seasonalised time series: https://docs.google.com/file/d/0B1Lm03a_91xia2RTOHZrajJtZXM/edit?usp=sharing

Below is the dput() of my data:

> dput(fr.monthly.temp.ts)
structure(c(2.7, 0.4, 4.7, 10, 13, 16.9, 19.2, 18.3, 15.7, 10.6,   
4.9, 3.5, 4.1, 3.2, 7.5, 10.3, 10, 15.1, 18.2, 17.4, 15, 10.2, 
6.3, 3.5, 3.8, 5.9, 7.6, 7.1, 12.9, 14.9, 17.6, 17.3, 15.5, 12.1, 
6.9, 2.7, 3, 4.6, 5.5, 10.3, 13.6, 16.3, 20.2, 18.5, 13.9, 11.2, 
5.4, 4.8, 1.7, 4, 7.4, 9.3, 11.9, 16.5, 20, 17.6, 14.7, 8.4, 
5.5, 3.8, 4.3, 3.1, 5.6, 8.5, 12.6, 16.1, 18.2, 18.9, 16, 12.7, 
7.4, 2.3, 2.5, 2.1, 6.3, 8.4, 12.7, 15.1, 16.5, 17.9, 16.2, 11.6, 
7.6, 5.6, 1.7, 4.8, 5, 7.7, 14.2, 16.8, 17.9, 17.1, 14.8, 12.1, 
6.5, 3.6, 2.2, 2, 4.7, 10.4, 12.8, 14.2, 16.3, 18, 14.2, 12.2, 
5, 4.9, 4, 5.4, 6.6, 8.5, 11.9, 16.1, 16.4, 17.3, 14.2, 11.9, 
5.9, 6, 1.6, 4.5, 6.4, 8.3, 13.6, 16.1, 20.8, 20.7, 17.5, 11.3, 
7.3, 6.6, 4.6, 6.8, 8.4, 9.2, 13.8, 15.5, 17.9, 15.5, 12.5, 10, 
5.5, 5.8, 5.4, 4.7, 7.9, 9.1, 13, 15.8, 16.5, 17.6, 15.4, 12.3, 
9.2, 4, 0.7, 6.5, 7.4, 11.2, 12.2, 15.3, 17.3, 18.2, 15.3, 10.6, 
6.3, 5.7, 3.5, 4.3, 5.7, 8.5, 14.2, 17, 17.2, 17.5, 14.7, 9.6, 
4.6, 7, 6.4, 4.8, 5.9, 9.5, 13.8, 14, 17.4, 18.4, 14.5, 11.5, 
7, 4.3, 1.1, 1.4, 4.4, 6.7, 15.1, 17.6, 18.3, 17.2, 16.4, 9.4, 
7.3, 1.4, 3.7, 5.4, 6.5, 8.4, 14.2, 15, 18, 18.1, 15.4, 9.7, 
6.4, 6.9, 3.3, 3.7, 6.2, 7.8, 13.8, 16.3, 15.9, 18.9, 16.2, 8.8, 
4.6, 5.5, 5, 6.4, 8.2, 9.9, 14.4, 16, 17.4, 16.5, 15.2, 11.5, 
6, 4, 6.4, 4.2, 7.2, 8.9, 13.7, 16.9, 20.6, 18, 17, 14.1, 4.7, 
4.5, 3.4, 4.7, 6.6, 8, 14.8, 16.3, 16.7, 16.9, 13.7, 9.2, 5.4, 
4.5, 3.7, 6.3, 7.6, 9.4, 12.2, 14.1, 19.9, 18.8, 15.1, 12.3, 
5.3, 3.8, 3.8, 2.4, 6.4, 9.2, 14.1, 16.2, 18, 15.9, 15.2, 11.7, 
7.1, 4.5, 4.8, 5.6, 4.3, 9.1, 12.9, 17, 18, 17.6, 13.3, 11.8, 
4.9, 3.9, 4.1, 8.3, 7.2, 10.3, 11.6, 14.5, 18.2, 18.7, 17.3, 
11.5, 8.3, 2.5, 4.3, 4.5, 7.7, 9.8, 13.7, 15.7, 18, 17.8, 15.2, 
11.3, 6.7, 2.9, 5, 6.4, 7.1, 9.3, 11.8, 16.1, 20.5, 19.3, 15.8, 
11.5, 8.2, 3.7, 0.3, -0.2, 6.7, 7.8, 13.2, 16.3, 19.1, 18.1, 
18.4, 11.4, 7.3, 6.4, 5.8, 3.3, 7, 9.7, 12.1, 17.7, 17.3, 18.2, 
15.9, 11.9, 8.6, 4.5, 3.7, 3.3, 5.8, 8.8, 13.8, 17.5, 17.7, 17, 
12.8, 10.6, 8.2, 3.2, 4.8, 1.4, 5.5, 8, 12.1, 15.8, 17.4, 20.4, 
17.2, 11, 7.4, 5, 1.8, 4.3, 7.8, 10.1, 13.1, 15.4, 19.5, 20.1, 
16.7, 12, 5.5, 0.3, 3.3, 3.1, 6.3, 10.4, 13.8, 17.2, 20, 17.5, 
17.1, 11.9, 5.8, 7.6, 2.6, 5.1, 6.2, 9.1, 11.6, 17.2, 19.5, 18.1, 
16.1, 10.7, 7, 3.9, 6.5, 4.6, 7.9, 8.3, 13.4, 16.1, 17.2, 18, 
16, 9.1, 6.6, 4.2, 5.3, 6.9, 5.6, 9.9, 14.2, 16.6, 18.6, 19.1, 
15.5, 11.7, 6.3, 3.2, 4.4, 3.9, 8.8, 7.7, 11.7, 16.8, 17.5, 18.2, 
15.6, 11.3, 9.3, 2.5, 5.3, 4.7, 5.4, 10.2, 11.5, 16.4, 17.3, 
18.1, 15.2, 10.3, 8.7, 2.6, -0.9, 4.5, 7.1, 9.6, 13.5, 17.1, 
17.1, 17.5, 15.6, 10.6, 7.6, 1.1, 0.7, 4.5, 7.3, 8.2, 10.3, 16.8, 
19.3, 16.9, 15.5, 10.8, 6.6, 3.7, -0.2, -0.1, 7.7, 10.6, 13.1, 
16.7, 18.1, 18.7, 16.7, 13.2, 5.5, 4.8, 4.8, 5.3, 8, 11.5, 14.2, 
16.4, 19.2, 19.2, 16, 12.4, 5.9, 3.4, 5.1, 2.2, 5.1, 11.1, 13.4, 
16, 18.6, 20.6, 15.2, 10.1, 7.1, 3.4, -1, 7.1, 8.4, 11.9, 14.8, 
17.8, 20, 18.1, 16.7, 12.3, 6.5, 4.8, 1.7, 6.4, 6.7, 11.2, 13.1, 
15.7, 18.9, 17.9, 16.2, 11.3, 7.1, 2.1, 1, 1.3, 7.3, 11.3, 14.8, 
17.9, 20.4, 20.9, 17.6, 12.1, 8.3, 3.8, 5.7, 4.5, 9.5, 10.4, 
14, 15.8, 17, 17.8, 15.5, 11.4, 7.2, 4.6, 4.5, 5.4, 5.7, 11.7, 
12.2, 16.8, 20.6, 19.8, 18.6, 13.4, 6.4, 5.1, 3, 6.4, 8, 8.7, 
14.2, 18.3, 20.2, 18.6, 15.2, 11.4, 7.4, 1.1, 4.6, 4.7, 5.8, 
9.1, 11.8, 16.1, 18.7, 17.5, 16.5, 10.5, 8.7, 4.9, 2.7, 2.8, 
8.1, 11.2, 14.5, 17.9, 20.2, 18.9, 13.1, 10.9, 5.5, 3.5, 1.1, 
3, 7.5, 10.1, 14.8, 15.4, 18, 18.8, 16.2, 12.1, 7, 6.8, 1.7, 
2.3, 7.5, 8.6, 12.6, 16, 16.4, 16.9, 15.5, 12.4, 8, 6.2, 4.4, 
3.6, 4.6, 10.3, 12.5, 16.4, 19.1, 19.2, 15.7, 10.4, 6.7, 6.4, 
4.4, -1.8, 6.7, 8.1, 13.8, 14.4, 17.8, 16.4, 16.4, 10.6, 5.3, 
5.2, 3.1, 6.9, 9.8, 9.6, 11.5, 17, 18.5, 17.6, 15.1, 11.8, 6.8, 
3.6, 3.7, 6.2, 4.9, 7.9, 13.9, 15.6, 17.9, 18.4, 17.3, 11.4, 
6.7, 5.1, 3.4, 4.5, 8.6, 10.2, 13.8, 17, 20.3, 18.9, 17.2, 12.2, 
6.8, 5.7, 3.5, 5, 8, 9.6, 14.5, 17.6, 16.8, 17.3, 14.5, 11.1, 
8.4, 3.5, 3.6, 7.6, 8.3, 11.7, 12.5, 16.6, 17.7, 18, 18.5, 12.3, 
6.4, 4.5, 4.8, 3.7, 3.9, 9.1, 11.5, 15.8, 17.6, 18.6, 15.5, 11.9, 
5.4, 1.3, -1.6, -0.3, 6.5, 9.6, 12.2, 15.8, 18.5, 16.5, 15.2, 
11.5, 9.3, 1.3, 1.5, 5.2, 5.6, 9.6, 14.5, 16.8, 19.6, 18.2, 16.7, 
9.6, 7.2, 3.2, 3.6, 1.7, 6.6, 8.7, 12.7, 16.1, 16.7, 17.1, 13.7, 
12.2, 6.3, 5.7, 2.6, 7.9, 6.2, 10.5, 13.2, 17, 16.8, 17.2, 16.6, 
12.7, 5, 5.3, 3.5, 5.5, 7.7, 8.8, 12.5, 15.6, 19.8, 18.1, 15.3, 
13.2, 7.1, 3, 3.3, 4.3, 6.8, 9.9, 11.8, 15.9, 17.8, 17.2, 15.1, 
13.5, 6.8, 3, 4.8, 2.1, 6.2, 9.2, 13.2, 15, 19.1, 18.1, 15.9, 
13.1, 7.1, 1.4, 4.1, 4.3, 4.4, 7.6, 12.8, 17.6, 17.8, 18.3, 16.6, 
11.3, 8.7, 2.6, 3.1, 4.2, 3.8, 10.5, 13.7, 14.8, 19.7, 18.7, 
15.7, 12.3, 5.8, 4.9, 3.2, 5.5, 7.9, 8.9, 11.7, 14.3, 18, 17.1, 
13.3, 10.9, 7.3, 4.5, 3, 3.4, 6.1, 7.6, 13.5, 17, 18.1, 19.9, 
16.7, 10.6, 6.8, 3.7, 6.2, 5.5, 7.3, 9.4, 12.5, 15.9, 17.7, 18.6, 
14.5, 8.2, 7.4, 6.8, 6.4, 5.5, 5.3, 9, 12.1, 15.9, 19.1, 19.8, 
16.1, 10.4, 6.7, 3.1, 4.1, 4.8, 6, 8.9, 14, 18.8, 20.1, 19, 14.8, 
11.8, 6.6, 3.1, 3.7, 6.6, 8.3, 8.3, 12.1, 14.8, 17.8, 16.9, 14.7, 
12.9, 7, 5.3, 3.3, 4, 7.2, 7.8, 12.3, 15.2, 17.3, 17.2, 15.6, 
11.8, 6.7, 5.1, 1.3, 4, 6.6, 8.2, 12.3, 16.5, 18.5, 17.1, 15.7, 
12.4, 6.7, 5.7, 2.2, 6.3, 6.2, 8.4, 11.9, 15, 16.4, 18.6, 16.5, 
10.8, 5.8, 3.1, 3.3, 2.9, 9.2, 10, 12.6, 16, 17.5, 18.8, 16.2, 
11.2, 7.2, 3.8, 4.6, 5, 6.3, 9.3, 13.4, 17.4, 20.1, 18, 17.4, 
11.4, 8.3, 4.9, 5.5, 2.5, 7, 8.9, 11.5, 17.1, 22.2, 19.3, 16.3, 
11.9, 7.6, 4.5, 4.2, 3.5, 5.2, 9.6, 10.4, 15.8, 18.8, 18.4, 14.7, 
11.9, 9, 4.5, -1, 3.8, 5.2, 9.7, 12.5, 15.3, 19.4, 17.6, 17.3, 
12.3, 4.4, 5.6, 3.9, -0.6, 5.9, 6.9, 13.7, 16.9, 18.7, 17.6, 
14.9, 13.1, 7.9, 5, -0.8, 3.7, 4.8, 10.9, 11.4, 15, 18.6, 18.6, 
17.8, 12.4, 7.1, 5.2, 6.4, 4.9, 6.5, 10.1, 13.8, 16.2, 17.8, 
18.7, 15.7, 12.9, 6.3, 6, 4.2, 5.6, 9.3, 8.2, 15.3, 16.9, 20.2, 
19.5, 16.5, 13.2, 7, 5.6, 4.8, 8.8, 8.7, 8.9, 15.3, 16, 19.7, 
20.4, 15.9, 13.3, 7.2, 3.1, 3.9, 1.9, 9, 8.7, 11.7, 14.9, 19.6, 
20.7, 17.9, 10.9, 6.9, 3.6, 2.8, 4.9, 7.6, 9.5, 15.3, 16.1, 19.1, 
19.9, 15.5, 9.6, 9, 4.8, 5.9, 3.5, 7, 10.4, 14.1, 17.3, 17.8, 
18.7, 14.7, 10.4, 4.8, 6.2, 5.2, 5.1, 9.4, 8.7, 13.6, 17.1, 21.4, 
19.9, 15, 12, 10.2, 6.5, 4.5, 7.5, 6.5, 9.9, 13.6, 16.1, 21.1, 
20.2, 14.5, 14.6, 7.5, 3.8, 5, 2.9, 6, 10, 12.2, 17.5, 18.7, 
18.2, 14.2, 11.9, 6.9, 3.4, 2.3, 6.9, 9.3, 10, 14.2, 16.3, 18.6, 
21, 17, 12.4, 8.4, 5.5, 5, 5.9, 8.1, 9, 14.9, 17, 18.5, 19.4, 
16.1, 11.6, 5.2, 4.5, 5.3, 4.3, 8, 10, 15.2, 16.3, 20.2, 19.4, 
17.9, 12.2, 6.4, 5, 3.7, 6.6, 7.5, 9.9, 15, 17.8, 17.5, 19.6, 
16.9, 12.2, 8.2, 7.1), .Tsp = c(1901, 2000.91666666667, 12), class = "ts")  

I run stl() on it to remove the seasonality:

# calculate and remove the seasonality  
fr.monthly.temp.ts.stl <- stl(fr.monthly.temp.ts, s.window="periodic")    # get the    components  
fr.monthly.temp.seas <- fr.monthly.temp.ts.stl$time.series[,"seasonal"]  
#plot(fr.monthly.temp.seas)  

fr.monthly.temp.ts.noseas <- fr.monthly.temp.ts - fr.monthly.temp.seas  
#plot(fr.monthly.temp.ts.noseas)  

Then remove the trend with a regression:

fr.mtrend.noseas <- lm(fr.monthly.temp.ts.noseas~t)  
summary(fr.mtrend.noseas)  

and then use the residuals of this model to fit an ARIMA model (after checking the ACF and PACF for which one is appropriate):

# create time series of residuals..this is our "detrended" series..for now use only linear trend result  
fr.monthly.temp.ts.new <- ts(fr.mtrend.noseas$resid, start=c(1901,1), frequency=12)
#plot.ts(fr.monthly.temp.ts.new, main="Detrended and de-seasonalized time series")

# ARIMA 1,1,1  
fit6 <- arima(fr.monthly.temp.ts.new,order=c(1,1,1))  
fit6  
tsdiag(fit6)  

I then make a prediction on the stationary time series:

#forecast for the stationary TS, for next 50 yrs months  
forecast <- predict(fit6,n.ahead=600)  

And then add back the trend and seasonality:

t.new <- (n+1):(n+600)  

#initial time series = stationaryTS + seasonality + trend  
fr.monthly.temp.ts.init <- fr.monthly.temp.ts.new + fr.monthly.temp.seas +
                            fr.mtrend.noseas$coefficients[1] + t * fr.mtrend.noseas$coefficients[2]  

#same for the prediction: we need to add seasonality and trend  
pred.Xt <- forecast$pred + fr.monthly.temp.seas[1:(1+50*12 - 1)] + 
                                fr.mtrend.noseas$coefficients[1] + t.new * fr.mtrend.noseas$coefficients[2]  

plot(fr.monthly.temp.ts.init,type="l",xlim=c(1940,2060))  
lines(pred.Xt,col="red",lwd=2)  

So going back to my question: Do I need to add some white noise to the prediction to be able to realistically predict temperature? And more generally, is my method correct?

$\endgroup$
1
$\begingroup$

I think you're doing this correctly. The predictions show the deterministic component of the model, as intended. As a loose analogy, if you add a trend line to a simple scatter plot (abline(lm(...), e.g.), you wouldn't expect the trend line to wobble around. Similarly, the forecast represents the best guess at the future temperatures.

Were you perhaps interested in a stochastic simulation, or bounding the forecast estimates with confidence intervals? Or maybe there is another seasonal component in the observed time series that is missing from your forecast? I'll elaborate on these three possibilities below.

  1. For the confidence intervals, look at forecast$se. Multiply that by the appropriate critical value, e.g., qnorm(0.975), and plot using lines(), and you can add some confidence bands.
  2. If you want a stochastic simulation, I would simulate the time series manually using the equation from the model that you fitted. Basically, simulate that prediction 1 time step at a time in a for() loop, and each time step add a Wiener noise term. For an AR(1) process with a coefficient of 1.00, this could amount to overlaying Brownian motion (diffinv(rnorm(40))). Your model has more terms, so it's not so simple. Try something akin to X[i] <- B1*X[i-1] + Sigma*rnorm(1), where Sigma is the s.d. of the noise that you want to add, B1 is the AR(1) coefficient, and X is the response variable. Adjust as necessary to include your seasonal and MA terms, etc. You could choose Sigma based on the residual variance of your fitted time series model.
  3. Perhaps what "looks wrong" to you is absence of a seasonality term with a multi-year period. For example, the effect of the North Atlantic Oscillation on temperature. You could try adding the NAO Index as a covariate to the model.

To directly answer your questions:

Do you need to add a white noise term? No. If it's noise you want to add, do a stochastic simulation, which would allow the noise to propagate naturally (a shock at time t-1 will have an effect at time t b/c X[t] is correlated with X[t-1]; also, the MA term means that the shocks/errors are correlated b/c Epsilon[t] is correlated with Epsilon[t-1]).

Is you method correct? Given the structure of your model, your forecast seems reasonable. If the forecast is lacking some deterministic pattern, try accounting for things like NAO. If you desire an explicit graphical representation of uncertainty (due to lack of additional terms, or other process/ observation errors), I would suggest communicating that the model upon which your forecast is based has some residual variance --- put s.e.'s or CI's on the forecast.

My hunch is that you're mostly interested in doing a stochastic simulation. That would be the way to add randomness.

$\endgroup$
  • $\begingroup$ Thanks for the answer. My comments: 1. I am in fact interested in the stochastic simulation that you propose. I did not plot them in my chart here but I do calculate the confidence bands and they follow the general shape of the prediction. So the for()loop you suggest in #2 looks like a good idea. However, would a SARIMA model perhaps capture the seasonal randomness, I remember reading it in my notes...? 2. I was also looking for a way to add the multi-yr seasonality since the ACF and PACF do show a significant spike at around 22 months (see edit). But i did not know how to do that.. $\endgroup$ – avg May 5 '13 at 8:20
  • $\begingroup$ arima(..., seasonal=list(order=c(P,D,Q), period=22) That should allow you to incorporate a 22-month seasonal trend. Also, you and I should make sure we're on the same page with terminology -- I don't consider the "seasonal pattern" to be "seasonal noise". If I can look at a time series and immediately notice deviations from the mean that are regularly repeating, I consider that noise, not pattern. If I see one deviation that is consistent with a known driver (e.g., time series of boat traffic shows a pulse on the date when a lake was stocked with trophy fish), that's pattern too. $\endgroup$ – rbatt May 5 '13 at 11:58
  • 1
    $\begingroup$ just as an update, I implemented the stochastic simulation today for a simpler AR(1) model and got the result that I wanted. I will now try to implement it in ARMA..thanks again for the answer.. $\endgroup$ – avg May 10 '13 at 17:42
  • $\begingroup$ No problem. I'm glad your closing in on the desired analysis! $\endgroup$ – rbatt May 10 '13 at 17:58
0
$\begingroup$

Your approach is to assume certain structure e.g. a linear trend and fixed monthly effects. The correct approach is to incorporate SARIMA and any level shifts , time trends , monthly effects that may be needed. Furtermore one needs to test for constancy of the variance of the error process and constancy of model parameters over time. The need for "randomness" in the prediction can be met by drawing errors from the set of errors and using them to peturb the expected value thus creating a simulation (one of many ! )

$\endgroup$
  • $\begingroup$ Thanks for the answer. I have implemented the suggestion by @rbatt to introduce the random features and I will now explore how to do the same using a SARIMA approach. $\endgroup$ – avg May 10 '13 at 17:40

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.