-
Notifications
You must be signed in to change notification settings - Fork 8
/
Copy pathlec2.Rmd
161 lines (139 loc) · 5 KB
/
lec2.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
## 5. Fitting and choosing models
### 5.1 Back to the temperature series
```{r}
load("meteo.RData")
# plot data and trend:
plot(T.outside~date,meteo,type='l')
meteo$T.per = 18.2-4.9*sin(pi*(meteo$hours+1.6)/12)
lines(T.per~date,meteo,col='red')
```
### 5.2 Removing the diurnal periodicity
Assuming this is a sinus function,
$$\alpha_1 + \alpha_2 \sin(t + \alpha_3)$$,
we need non-linear regression ($\alpha_3$)
```{r}
attach(meteo) # so we can use variable names directly:
f = function(x) sum((T.outside - (x[1]+x[2]*sin(pi*(hours+x[3])/12)))^2)
nlm(f,c(0,0,0))
```
### 5.3 Temperature anomaly
```{r}
plot(T.outside-T.per~date, meteo, type='l')
title("anomaly")
```
### 5.4 What can we do with such models?
* try to find out which model fits best (model selection)
* learn how they were/could have been generated
* predict future observations (estimation/prediction/forecasting)
* generate similar data ourselves (simulation)
### 5.5 How to select a ``best'' model?
A possible approach is to find the minimum for {\em Akaike's Information
Criterion} (AIC) for ARMA($p,q$) models and series of length $n$:
$$AIC = \log \hat{\sigma}^2 + 2(p+q+1)/n$$
with $\hat{\sigma}^2$ the estimated residual (noise) variance.
Instead of finding a single best model using this single criterion, it
may be better is to select a small group of ``best'' models, and look
at model diagnostics for each: is the residual white noise? does it have
stationary variance?
Even better may be to keep a number of ``fit'' models and consider each as
(equally?) suitable candidates.
### 5.6 AIC for AR(p), and as as a function of $p$
```{r}
n = 10
aic = numeric(n)
for (i in 1:n)
aic[i] = arima(T.outside, c(i,0,0))$aic # AR(i)
aic
plot(aic[2:10], type='l')
```
### 5.7 Anomaly AIC as a function of $p$, for AR(p)
```{r}
an = T.outside - T.per
aic_an = numeric(n)
for (i in 1:n)
aic_an[i] = arima(an,c(i,0,0))$aic # AR(i)
aic_an
plot(aic_an[2:10], type='l')
```
```{r}
# Prediction, e.g. with AR(6)
x = arima(T.outside, c(6,0,0))
pltpred = function(xlim, Title) {
plot(T.outside, xlim = xlim, type='l')
start = nrow(meteo) + 1
pr = predict(x, n.ahead = xlim[2] - start + 1)
lines(start:xlim[2], pr$pred, col='red')
lines(start:xlim[2], pr$pred+2*pr$se, col='green')
lines(start:xlim[2], pr$pred-2*pr$se, col='green')
title(Title)
}
pltpred(c(9860, 9900), "10 minutes")
pltpred(c(9400, 10000), "110 minutes")
pltpred(c(8000, 11330), "predicting 1 day")
pltpred(c(1, 19970), "predicting 1 week")
plot(T.outside,xlim=c(1,19970), type='l')
x.an = arima(an, c(6,0,0)) # model the anomaly by AR(6)
x.pr = as.numeric(predict(x.an, 10080)$pred)
x.se = as.numeric(predict(x.an, 10080)$se)
hours.all = c(meteo$hours, max(meteo$hours) + (1:10080)/60)
T.per = 18.2-4.9*sin(pi*(hours.all+1.6)/12)
lines(T.per, col = 'blue')
hours.pr = c(max(meteo$hours) + (1:10080)/60)
T.pr = 18.2-4.9*sin(pi*(hours.pr+1.6)/12)
lines(9891:19970, T.pr+x.pr, col='red')
lines(9891:19970, T.pr+x.pr+2*x.se, col='green')
lines(9891:19970, T.pr+x.pr-2*x.se, col='green')
title("predicting 1 week: periodic trend + ar(6) for anomaly")
x = arima(T.outside, c(6,0,0))
plot(T.pr + arima.sim(list(ar = x.an$coef[1:6]), 10080, sd = sqrt(0.002556)),
col = 'red', ylab="Temperature")
lines(18+arima.sim(list(ar = x$coef[1:6]), 10080, sd=sqrt(0.002589)),
col = 'blue')
title("red: with trend, blue: without trend")
```
### 5.8 What can we learn from this?
Prediction/forecasting:
* AR(6) prediction is a compromise between the end of the series and the trend
* the closer we are to observations, the more similar the prediction is to the
nearest (last) observation
* further in the future the prediction converges to the trend
* the more useful (realistic) the trend is, the more realistic
the far-into-the-future prediction becomes
* the standard error of prediction increases when predictions are further in the
future.
### 5.9 A side note: fitting a phase with a linear model
A phase shift model $\alpha \sin(x + \phi)$ can also be modelled by
$\alpha sin(x) + \beta cos(x)$; this is essentially a re-parameterization.
Try the following code:
```{r eval=FALSE}
x = seq(0, 4*pi, length.out = 200) # x plotting range
f = function(dir, x) { # plot the combination of a sin+cos function, based on dir
a = sin(dir)
b = cos(dir)
# three curves:
plot(a * sin(x) + b * cos(x) ~ x, type = 'l', asp=1, col = 'green')
lines(x, a * sin(x), col = 'red')
lines(x, b * cos(x), col = 'blue')
# legend:
lines(c(10, 10+a), c(2,2), col = 'red')
lines(c(10, 10), c(2,2+b), col = 'blue')
arrows(x0 = 10, x1 = 10+a, y0 = 2, y1 = 2+b, .1, col = 'green')
title("red: sin(x), blue: cos(x), green: sum")
}
for (d in seq(0, 2*pi, length.out = 100)) {
f(d, x)
Sys.sleep(.1) # give us some time to think this over!
}
```
So, we can fit the same model by a different parameterization:
```{r}
nlm(f,c(0,0,0))$minimum
tt = pi * hours / 12
g = function(x) sum((T.outside - (x[1]+x[2]*sin(tt)+x[3]*cos(tt)))^2)
nlm(g,c(0,0,0))
```
which is a linear model, identical to:
```{r}
lm(T.outside~sin(tt)+cos(tt))
```
### 5.10 Next steps: how do we fit coefficients?