Predictions of a monthly temperature time series: adding noise to the predicted values

Question

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?

Answer 1

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.

Answer 2

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 ! )