Practical Time Series Analysis
Class notes
Practical Time Series Analysis
- Introducción
- Semana 1 - Revisión de estadisticos básicos e inferencia
- Semana 2 - Visualizando series temporales, y comenzando a modelas series temporales
- Semana 3 - Estacionariedad, procesos MA(q) y AR (p)
- Semana 4 - procesos AR(p), ecuaciones Yule-Walker y PACF
- Semana 5. Akaike Information Criterion (AIC), Mixed Models, Integrated Models
- Semana 6 - Estacionalidad, SARIMA y Predicción (Forecasting)
Introducción
Instalación de R y RStudio. Instalar paquete faraway, para los datasets, y fpp2 para el manejo de series temporales.
Primer ejercicio
- Carga el dataset de venenos, tratamientos y ratas (efecto del agente tóxico en ratas).
data("rats")- Escribe un párrafo explicando las variables que tenemos en el dataset y de donde Faraway cogio los datos.
glimpse(rats)## Observations: 48
## Variables: 3
## $ time <dbl> 0.31, 0.82, 0.43, 0.45, 0.45, 1.10, 0.45, 0.71, 0.46, 0...
## $ poison <fct> I, I, I, I, I, I, I, I, I, I, I, I, I, I, I, I, II, II,...
## $ treat <fct> A, B, C, D, A, B, C, D, A, B, C, D, A, B, C, D, A, B, C...summary(rats)## time poison treat
## Min. :0.1800 I :16 A:12
## 1st Qu.:0.3000 II :16 B:12
## Median :0.4000 III:16 C:12
## Mean :0.4794 D:12
## 3rd Qu.:0.6225
## Max. :1.2400Tenemos el timestamp como double, la variable poison como un factor, y el tratamiento como otro factor.
- Obten un box plot para los tiempos de supervivencia por veneno y por tratamiento.
plot(time~treat,data=rats) +
title("Suvival time per treatment")
## integer(0)plot(time~poison,data=rats) +
title("Survival time per poison")
## integer(0)Semana 1 - Revisión de estadisticos básicos e inferencia
Concatenación, resumen de 5 numeros y desviación estándar
data.1=c(35, 8, 10, 23, 42)
data.1## [1] 35 8 10 23 42print(data.1)## [1] 35 8 10 23 42summary(data.1)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 8.0 10.0 23.0 23.6 35.0 42.0mean(data.1)## [1] 23.6sd(data.1)## [1] 14.97665sum(data.1)/5## [1] 23.6sum(data.1)/length(data.1)## [1] 23.6Histograma
small.size.dataset=c(91,49,76,112,97,42,70, 100, 8, 112, 95, 90, 78, 62, 56, 94, 65, 58, 109, 70, 109, 91, 71, 76, 68, 62, 134, 57, 83, 66)
hist(small.size.dataset)
hist(small.size.dataset, xlab='My data points')
hist(small.size.dataset, xlab='My data points', main='Histogram of my data')
hist(small.size.dataset, xlab='My data points', main='Histogram of my data', freq=F)
hist(small.size.dataset, xlab='My data points', main='Histogram of my data', freq=F, col='green')
hist(small.size.dataset, xlab='My data points', main='Histogram of my data', freq=F, col='green')
lines(density(small.size.dataset))
hist(small.size.dataset, xlab='My data points', main='Histogram of my data', freq=F, col='green')
lines(density(small.size.dataset), col='red', lwd=5)
hist(small.size.dataset, xlab='My data points', main='Histogram of my data', freq=F, col='green', breaks=10)
hist(small.size.dataset, xlab='My data points', main='Histogram of my data', freq=F, col='green', breaks=10)
lines(density(small.size.dataset), col='red', lwd=5)
ScatterPlot
set.seed(2016) # There is a typo in the video (set.seed=2016)
Test_1_scores=round(rnorm(50, 78, 10))
Test_2_scores=round(rnorm(50, 78, 14))
Test_1_scores # Data won't be the same with the data generated in the video lecture since there was a typo in set.seed. ## [1] 69 88 77 81 50 75 70 71 82 80 74 75 91 69 94 81 69
## [18] 88 71 71 74 71 75 84 80 68 72 91 89 67 68 77 90 63
## [35] 79 92 70 70 73 63 77 83 75 72 87 75 65 100 86 61Test_2_scores # Data won't be the same with the data generated in the video lecture since there was a typo in set.seed. ## [1] 76 66 43 62 80 55 101 61 73 72 88 87 70 70 85 97 74
## [18] 62 72 90 91 83 51 78 77 89 81 91 78 73 86 85 86 67
## [35] 57 106 70 87 74 81 97 86 61 88 73 70 80 74 46 95plot(Test_2_scores~Test_1_scores)
plot(Test_2_scores~Test_1_scores, main='Test scores for two exams (50 students)', xlab='Test_1_scores', ylab='Test 2 scores')
plot(Test_2_scores~Test_1_scores, main='Test scores for two exams (50 students)', xlab='Test_1_scores', ylab='Test 2 scores', col='blue')
Basic statistics I y II - Simple lineal regression
Vamos a usar la serie de datos de CO2 incluida en R para mostrar como hacer una regresión lineal y medir su ajuste.
Primero damos un vistazo rápido a la serie y la dibujamos.
print(co2)## Jan Feb Mar Apr May Jun Jul Aug Sep Oct
## 1959 315.42 316.31 316.50 317.56 318.13 318.00 316.39 314.65 313.68 313.18
## 1960 316.27 316.81 317.42 318.87 319.87 319.43 318.01 315.74 314.00 313.68
## 1961 316.73 317.54 318.38 319.31 320.42 319.61 318.42 316.63 314.83 315.16
## 1962 317.78 318.40 319.53 320.42 320.85 320.45 319.45 317.25 316.11 315.27
## 1963 318.58 318.92 319.70 321.22 322.08 321.31 319.58 317.61 316.05 315.83
## 1964 319.41 320.07 320.74 321.40 322.06 321.73 320.27 318.54 316.54 316.71
## 1965 319.27 320.28 320.73 321.97 322.00 321.71 321.05 318.71 317.66 317.14
## 1966 320.46 321.43 322.23 323.54 323.91 323.59 322.24 320.20 318.48 317.94
## 1967 322.17 322.34 322.88 324.25 324.83 323.93 322.38 320.76 319.10 319.24
## 1968 322.40 322.99 323.73 324.86 325.40 325.20 323.98 321.95 320.18 320.09
## 1969 323.83 324.26 325.47 326.50 327.21 326.54 325.72 323.50 322.22 321.62
## 1970 324.89 325.82 326.77 327.97 327.91 327.50 326.18 324.53 322.93 322.90
## 1971 326.01 326.51 327.01 327.62 328.76 328.40 327.20 325.27 323.20 323.40
## 1972 326.60 327.47 327.58 329.56 329.90 328.92 327.88 326.16 324.68 325.04
## 1973 328.37 329.40 330.14 331.33 332.31 331.90 330.70 329.15 327.35 327.02
## 1974 329.18 330.55 331.32 332.48 332.92 332.08 331.01 329.23 327.27 327.21
## 1975 330.23 331.25 331.87 333.14 333.80 333.43 331.73 329.90 328.40 328.17
## 1976 331.58 332.39 333.33 334.41 334.71 334.17 332.89 330.77 329.14 328.78
## 1977 332.75 333.24 334.53 335.90 336.57 336.10 334.76 332.59 331.42 330.98
## 1978 334.80 335.22 336.47 337.59 337.84 337.72 336.37 334.51 332.60 332.38
## 1979 336.05 336.59 337.79 338.71 339.30 339.12 337.56 335.92 333.75 333.70
## 1980 337.84 338.19 339.91 340.60 341.29 341.00 339.39 337.43 335.72 335.84
## 1981 339.06 340.30 341.21 342.33 342.74 342.08 340.32 338.26 336.52 336.68
## 1982 340.57 341.44 342.53 343.39 343.96 343.18 341.88 339.65 337.81 337.69
## 1983 341.20 342.35 342.93 344.77 345.58 345.14 343.81 342.21 339.69 339.82
## 1984 343.52 344.33 345.11 346.88 347.25 346.62 345.22 343.11 340.90 341.18
## 1985 344.79 345.82 347.25 348.17 348.74 348.07 346.38 344.51 342.92 342.62
## 1986 346.11 346.78 347.68 349.37 350.03 349.37 347.76 345.73 344.68 343.99
## 1987 347.84 348.29 349.23 350.80 351.66 351.07 349.33 347.92 346.27 346.18
## 1988 350.25 351.54 352.05 353.41 354.04 353.62 352.22 350.27 348.55 348.72
## 1989 352.60 352.92 353.53 355.26 355.52 354.97 353.75 351.52 349.64 349.83
## 1990 353.50 354.55 355.23 356.04 357.00 356.07 354.67 352.76 350.82 351.04
## 1991 354.59 355.63 357.03 358.48 359.22 358.12 356.06 353.92 352.05 352.11
## 1992 355.88 356.63 357.72 359.07 359.58 359.17 356.94 354.92 352.94 353.23
## 1993 356.63 357.10 358.32 359.41 360.23 359.55 357.53 355.48 353.67 353.95
## 1994 358.34 358.89 359.95 361.25 361.67 360.94 359.55 357.49 355.84 356.00
## 1995 359.98 361.03 361.66 363.48 363.82 363.30 361.94 359.50 358.11 357.80
## 1996 362.09 363.29 364.06 364.76 365.45 365.01 363.70 361.54 359.51 359.65
## 1997 363.23 364.06 364.61 366.40 366.84 365.68 364.52 362.57 360.24 360.83
## Nov Dec
## 1959 314.66 315.43
## 1960 314.84 316.03
## 1961 315.94 316.85
## 1962 316.53 317.53
## 1963 316.91 318.20
## 1964 317.53 318.55
## 1965 318.70 319.25
## 1966 319.63 320.87
## 1967 320.56 321.80
## 1968 321.16 322.74
## 1969 322.69 323.95
## 1970 323.85 324.96
## 1971 324.63 325.85
## 1972 326.34 327.39
## 1973 327.99 328.48
## 1974 328.29 329.41
## 1975 329.32 330.59
## 1976 330.14 331.52
## 1977 332.24 333.68
## 1978 333.75 334.78
## 1979 335.12 336.56
## 1980 336.93 338.04
## 1981 338.19 339.44
## 1982 339.09 340.32
## 1983 340.98 342.82
## 1984 342.80 344.04
## 1985 344.06 345.38
## 1986 345.48 346.72
## 1987 347.64 348.78
## 1988 349.91 351.18
## 1989 351.14 352.37
## 1990 352.69 354.07
## 1991 353.64 354.89
## 1992 354.09 355.33
## 1993 355.30 356.78
## 1994 357.59 359.05
## 1995 359.61 360.74
## 1996 360.80 362.38
## 1997 362.49 364.34help(co2)
class(co2)## [1] "ts"#dibujando el grafico
plot(co2,main="Concentración CO2 atmosférica")
Un modelo de regresión lineal sigue la expresion Yt=(B0 + B1*Xt) + Ruido. El ruido, que puede darse con otros nombres, suele tener media 0, una cierta varianza y el valor del ruido de un tiempo t no está correlado con el valor del ruido en otros momentos de tiempo. Es decir, se considera ruido blanco.
Creamos un modelo lineal y lo comparamos con la gráfica de los valores reales.
co2.linear.model = lm(co2 ~ time(co2))
co2.linear.model##
## Call:
## lm(formula = co2 ~ time(co2))
##
## Coefficients:
## (Intercept) time(co2)
## -2249.774 1.307#dibujamos la grafica original con la recta de regresión
plot(co2,main="Concentración CO2 atmosférica con recta de regresión")
abline(co2.linear.model)
#calculamos los residuos
co2.residuals <- resid(co2.linear.model)
plot(co2.residuals)
#comprobamos que los residuos son una distribución normal con media cero
hist(co2.residuals)
#visualizamos varios gráficos juntos
par(mfrow=c(1,3))
( c02.residuals = resid( co2.linear.model ) )## 1 2 3 4 5
## 3.808180799 4.589222742 4.670264684 5.621306627 6.082348569
## 6 7 8 9 10
## 5.843390511 4.124432454 2.275474396 1.196516339 0.587558281
## 11 12 13 14 15
## 1.958600224 2.619642166 3.350684109 3.781726051 4.282767993
## 16 17 18 19 20
## 5.623809936 6.514851878 5.965893821 4.436935763 2.057977706
## 21 22 23 24 25
## 0.209019648 -0.219938409 0.831103533 1.912145475 2.503187418
## 26 27 28 29 30
## 3.204229360 3.935271303 4.756313245 5.757355188 4.838397130
## 31 32 33 34 35
## 3.539439072 1.640481015 -0.268477043 -0.047435100 0.623606842
## 36 37 38 39 40
## 1.424648785 2.245690727 2.756732670 3.777774612 4.558816554
## 41 42 43 44 45
## 4.879858497 4.370900439 3.261942382 0.952984324 -0.295973733
## 46 47 48 49 50
## -1.244931791 -0.093889848 0.797152094 1.738194036 1.969235979
## 51 52 53 54 55
## 2.640277921 4.051319864 4.802361806 3.923403749 2.084445691
## 56 57 58 59 60
## 0.005487633 -1.663470424 -1.992428482 -1.021386539 0.159655403
## 61 62 63 64 65
## 1.260697346 1.811739288 2.372781231 2.923823173 3.474865115
## 66 67 68 69 70
## 3.035907058 1.466949000 -0.372009057 -2.480967115 -2.419925172
## 71 72 73 74 75
## -1.708883230 -0.797841287 -0.186799345 0.714242597 1.055284540
## 76 77 78 79 80
## 2.186326482 2.107368425 1.708410367 0.939452310 -1.509505748
## 81 82 83 84 85
## -2.668463806 -3.297421863 -1.846379921 -1.405337978 -0.304296036
## 86 87 88 89 90
## 0.556745907 1.247787849 2.448829792 2.709871734 2.280913676
## 91 92 93 94 95
## 0.821955619 -1.327002439 -3.155960496 -3.804918554 -2.223876611
## 96 97 98 99 100
## -1.092834669 0.098207274 0.159249216 0.590291158 1.851333101
## 101 102 103 104 105
## 2.322375043 1.313416986 -0.345541072 -2.074499129 -3.843457187
## 106 107 108 109 110
## -3.812415245 -2.601373302 -1.470331360 -0.979289417 -0.498247475
## 111 112 113 114 115
## 0.132794468 1.153836410 1.584878353 1.275920295 -0.053037763
## 116 117 118 119 120
## -2.191995820 -4.070953878 -4.269911935 -3.308869993 -1.837828050
## 121 122 123 124 125
## -0.856786108 -0.535744165 0.565297777 1.486339719 2.087381662
## 126 127 128 129 130
## 1.308423604 0.379465547 -1.949492511 -3.338450568 -4.047408626
## 131 132 133 134 135
## -3.086366684 -1.935324741 -1.104282799 -0.283240856 0.557801086
## 136 137 138 139 140
## 1.648843029 1.479884971 0.960926914 -0.468031144 -2.226989202
## 141 142 143 144 145
## -3.935947259 -4.074905317 -3.233863374 -2.232821432 -1.291779489
## 146 147 148 149 150
## -0.900737547 -0.509695605 -0.008653662 1.022388280 0.553430223
## 151 152 153 154 155
## -0.755527835 -2.794485892 -4.973443950 -4.882402007 -3.761360065
## 156 157 158 159 160
## -2.650318123 -2.009276180 -1.248234238 -1.247192295 0.623849647
## 161 162 163 164 165
## 0.854891590 -0.234066468 -1.383024525 -3.211982583 -4.800940641
## 166 167 168 169 170
## -4.549898698 -3.358856756 -2.417814813 -1.546772871 -0.625730928
## 171 172 173 174 175
## 0.005311014 1.086352956 1.957394899 1.438436841 0.129478784
## 176 177 178 179 180
## -1.529479274 -3.438437331 -3.877395389 -3.016353446 -2.635311504
## 181 182 183 184 185
## -2.044269562 -0.783227619 -0.122185677 0.928856266 1.259898208
## 186 187 188 189 190
## 0.310940151 -0.868017907 -2.756975964 -4.825934022 -4.994892080
## 191 192 193 194 195
## -4.023850137 -3.012808195 -2.301766252 -1.390724310 -0.879682367
## 196 197 198 199 200
## 0.281359575 0.832401517 0.353443460 -1.455514598 -3.394472655
## 201 202 203 204 205
## -5.003430713 -5.342388770 -4.301346828 -3.140304885 -2.259262943
## 206 207 208 209 210
## -1.558221001 -0.727179058 0.243862884 0.434904827 -0.214053231
## 211 212 213 214 215
## -1.603011288 -3.831969346 -5.570927403 -6.039885461 -4.788843519
## 216 217 218 219 220
## -3.517801576 -2.396759634 -2.015717691 -0.834675749 0.426366194
## 221 222 223 224 225
## 0.987408136 0.408450078 -1.040507979 -3.319466037 -4.598424094
## 226 227 228 229 230
## -5.147382152 -3.996340209 -2.665298267 -1.654256324 -1.343214382
## 231 232 233 234 235
## -0.202172440 0.808869503 0.949911445 0.720953388 -0.738004670
## 236 237 238 239 240
## -2.706962727 -4.725920785 -5.054878842 -3.793836900 -2.872794958
## 241 242 243 244 245
## -1.711753015 -1.280711073 -0.189669130 0.621372812 1.102414755
## 246 247 248 249 250
## 0.813456697 -0.855501361 -2.604459418 -4.883417476 -5.042375533
## 251 252 253 254 255
## -3.731333591 -2.400291648 -1.229249706 -0.988207763 0.622834179
## 256 257 258 259 260
## 1.203876121 1.784918064 1.385960006 -0.332998051 -2.401956109
## 261 262 263 264 265
## -4.220914166 -4.209872224 -3.228830281 -2.227788339 -1.316746397
## 266 267 268 269 270
## -0.185704454 0.615337488 1.626379431 1.927421373 1.158463316
## 271 272 273 274 275
## -0.710494742 -2.879452800 -4.728410857 -4.677368915 -3.276326972
## 276 277 278 279 280
## -2.135285030 -1.114243087 -0.353201145 0.627840798 1.378882740
## 281 282 283 284 285
## 1.839924682 0.950966625 -0.457991433 -2.796949490 -4.745907548
## 286 287 288 289 290
## -4.974865605 -3.683823663 -2.562781721 -1.791739778 -0.750697836
## 291 292 293 294 295
## -0.279655893 1.451386049 2.152427992 1.603469934 0.164511877
## 296 297 298 299 300
## -1.544446181 -4.173404239 -4.152362296 -3.101320354 -1.370278411
## 301 302 303 304 305
## -0.779236469 -0.078194526 0.592847416 2.253889359 2.514931301
## 306 307 308 309 310
## 1.775973243 0.267015186 -1.951942872 -4.270900929 -4.099858987
## 311 312 313 314 315
## -2.588817044 -1.457775102 -0.816733160 0.104308783 1.425350725
## 316 317 318 319 320
## 2.236392668 2.697434610 1.918476553 0.119518495 -1.859439562
## 321 322 323 324 325
## -3.558397620 -3.967355678 -2.636313735 -1.425271793 -0.804229850
## 326 327 328 329 330
## -0.243187908 0.547854035 2.128895977 2.679937920 1.910979862
## 331 332 333 334 335
## 0.192021804 -1.946936253 -3.105894311 -3.904852368 -2.523810426
## 336 337 338 339 340
## -1.392768483 -0.381726541 -0.040684599 0.790357344 2.251399286
## 341 342 343 344 345
## 3.002441229 2.303483171 0.454525114 -1.064432944 -2.823391001
## 346 347 348 349 350
## -3.022349059 -1.671307117 -0.640265174 0.720776768 1.901818711
## 351 352 353 354 355
## 2.302860653 3.553902596 4.074944538 3.545986481 2.037028423
## 356 357 358 359 360
## -0.021929635 -1.850887692 -1.789845750 -0.708803807 0.452238135
## 361 362 363 364 365
## 1.763280078 1.974322020 2.475363962 4.096405905 4.247447847
## 366 367 368 369 370
## 3.588489790 2.259531732 -0.079426325 -2.068384383 -1.987342440
## 371 372 373 374 375
## -0.786300498 0.334741444 1.355783387 2.296825329 2.867867272
## 376 377 378 379 380
## 3.568909214 4.419951157 3.380993099 1.872035042 -0.146923016
## 381 382 383 384 385
## -2.195881074 -2.084839131 -0.543797189 0.727244754 1.138286696
## 386 387 388 389 390
## 2.069328639 3.360370581 4.701412523 5.332454466 4.123496408
## 391 392 393 394 395
## 1.954538351 -0.294419707 -2.273377764 -2.322335822 -0.901293879
## 396 397 398 399 400
## 0.239748063 1.120790005 1.761831948 2.742873890 3.983915833
## 401 402 403 404 405
## 4.384957775 3.865999718 1.527041660 -0.601916397 -2.690874455
## 406 407 408 409 410
## -2.509832513 -1.758790570 -0.627748628 0.563293315 0.924335257
## 411 412 413 414 415
## 2.035377200 3.016419142 3.727461084 2.938503027 0.809544969
## 416 417 418 419 420
## -1.349413088 -3.268371146 -3.097329203 -1.856287261 -0.485245318
## 421 422 423 424 425
## 0.965796624 1.406838566 2.357880509 3.548922451 3.859964394
## 426 427 428 429 430
## 3.021006336 1.522048279 -0.646909779 -2.405867837 -2.354825894
## 431 432 433 434 435
## -0.873783952 0.477257991 1.298299933 2.239341876 2.760383818
## 436 437 438 439 440
## 4.471425761 4.702467703 4.073509645 2.604551588 0.055593530
## 441 442 443 444 445
## -1.443364527 -1.862322585 -0.161280642 0.859761300 2.100803243
## 446 447 448 449 450
## 3.191845185 3.852887127 4.443929070 5.024971012 4.476012955
## 451 452 453 454 455
## 3.057054897 0.788096840 -1.350861218 -1.319819276 -0.278777333
## 456 457 458 459 460
## 1.192264609 1.933306552 2.654348494 3.095390437 4.776432379
## 461 462 463 464 465
## 5.107474322 3.838516264 2.569558206 0.510600149 -1.928357909
## 466 467 468
## -1.447315966 0.103725976 1.844767919hist(co2.residuals, main= "Histogram of CO2 Residuals")
qqnorm(c02.residuals, main= "Normal Probability Plot") #si los residuos siguen una distribución normal, este gráfico debe dar una linea recta
qqline(c02.residuals)
plot(c02.residuals ~ time(co2), main="Residuals on Time")
par(mfrow=c(1,1))
plot(co2.residuals ~ time(co2), xlim = c(1960,1963),main="Zoom en los residuos")#vemos que la distribución de los residuos no es exactamente normal porque la componente oscilatoria se puede ver también en los residuos en este grafico de solo una parte pequeña, y esta oscilación quita aleatoriedad a los residuos y por eso la distribución no es exactamente normal
Basic statistics III - Inference
Usaremos el dataset sleep. Primero hacemos un gráfico con los datos usando un box plot.
plot(extra ~ group, data=sleep, main="Extra sleep in gosset data by group") #las medias están bastante separadas
También podemos evaluar como de relacionadas o no están dos variables con el test de t-Student. Vamos a aplicarlo al dataset sleep, que contiene los datos de tiempo de sueño de varias personas en respuesta a dos drogas distintas.
attach(sleep)
extra.1 <- extra[group==1]
extra.2 <- extra[group==2]
#ejecutamos el test t-Student
t.test(extra.1, extra.2, paired=TRUE, alternative="two.sided")##
## Paired t-test
##
## data: extra.1 and extra.2
## t = -4.0621, df = 9, p-value = 0.002833
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -2.4598858 -0.7001142
## sample estimates:
## mean of the differences
## -1.58El valor alto y negativo de t, asi como el valor muy cercano a cero de p (el valor aceptado como limite superior es 0.05) hace pensar que estas variables son diferentes estadisticamente, como se intuia en el box plot inicial de forma gráfica.
Basic statistics IV
Vamos a visualizar pares de variables para explorar posibles correlaciones entre ellas. También vamos a aprender a interpretar la correlación y la covarianza.
help(trees)
pairs(trees,pch=21,bg=c("red")) #visualizamos las variables por pares, para detectar cuales son las que mas correlación tienen
#calculemos la covarianza entre variables
cov(trees)## Girth Height Volume
## Girth 9.847914 10.38333 49.88812
## Height 10.383333 40.60000 62.66000
## Volume 49.888118 62.66000 270.20280#y tambien calculamos la correlación entre variables
cor(trees)## Girth Height Volume
## Girth 1.0000000 0.5192801 0.9671194
## Height 0.5192801 1.0000000 0.5982497
## Volume 0.9671194 0.5982497 1.0000000#vemos que tanto la covarianza entre girth y volume asi como la correlación nos indican que están relacionadasSemana 2 - Visualizando series temporales, y comenzando a modelas series temporales
Visualizando y describiendo series temporales: graficos de tiempo, autocovarianza y autocorrelación
Vamos a usar el paquete astsa
Package by Robert H. Shumway and David S. Stoffer Contains data sets and scripts to accompany “Time Series Analysis and Its Applications: With R Examples” Aqui en enlace a la documentación del paquete
Una serie temporal es una colección de datos recogidos a lo largo del tiempo. Recoger muestras adyacentes en el tiempo introduce una cierta correlación entre estas muestras contiguas en el tiempo.
Series incluidas en el paquete astsa:
- jj - serie temporal de beneficios trimestrales de Johnson & Johnson. Ya es una serie temporal, no necesitamos usar la función ts() para convertir el dataset en una serie temporal. Con el grafico podemos ver cierta tendencia, asi como algunas fluctuaciones (algunas estacionales). También las variaciones se han agrandado con el tiempo.
- flu - numero de muertes mensuales por gripe y neumonía cada 10.000 habitantes en Estados Unidos.
- globtemp - desviaciones de la temperatura media de la superficie de los oceanos.
- globtempl - desviaciones de la temperatura media de la superficie terrestre.
#grafico sencillo serie JJ
plot(jj,type='o',main="J&J ganacias trimestrales por acción",xlab="Años",ylab="Beneficios")
#grafico sencillo serie flu
plot(flu,main="Muertes mensuales por gripe y neumonia en USA",xlab="Meses",ylab="Numero muertes")
#grafico sencillo serie globtemp
plot(globtemp,type='o',main="Desviaciones de la media de la temperatura de la superficie de los oceanos",xlab="Años",ylab="Desviación de la temperatura")
#grafico sencillo serie globtempl
plot(globtempl,type='o',main="Desviaciones de la media de la temperatura de la superficie terrestre",xlab="Años",ylab="Desviación de la temperatura")
#grafico sencillo serie star
plot(star,type='o',main="Magnitud de una estrella durante 600 dias",xlab="Dias",ylab="Magnitud")
Estacionariedad
Una serie temporal es estacionaria si no hay cambios significativos en la media (eso es, no hay una tendencia) ni en la varianza ni hay fluctuaciones periódicas. Esto significa que las propiedades de una parte de la serie son iguales a las propiedades de cualquier otra parte de la serie.
Hay un tipo de estacionariedad débil, donde no hay cambios en la media, aunque si que hay cambios en la varianza y hay estacionalidad.
Para poder modelar una serie temporal debe ser estacionaria (al menos débil).
Para series no estacionales (que son practicamente todas las que describen un proceso real) debemos aplicar transformaciones para hacerlas estacionarias.
Función de autocovarianza
La función de covarianza entre dos variables mide la dependencia lineal que hay entre esas variables. La función de covarianza es commutativa.
Un proceso estocastico es una colección de variables estadísticas. Podemos definir una serie temporal como la realización de un proceso estocástico.
La función de autocovarianza en nuestra serie temporal (realización de un proceso estocástico) se define como la covarianza entre las diferentes realizaciones (o puntos de la serie temporal).
Si la serie temporal es estacionaria, podemos decir la la covarianza entre un termino Xt y X(t+k) sera siempre igual y que solo depende de la diferencia temporal entre ambas muestras, no del momento temporal de cada muestra. Esta función se llama ɣ(k) o también coeficientes de autocovarianza, y solo se puede aplicar siempre que la serie temporal sea estacionaria.
En R la función de autocovarianza se calcula con la función acf(serie_temporal, type=‘covariance’). vamos a probarlo con una serie temporal totalmente aleatoria.
aleatorio <- ts(rnorm(100)) #genera una serie temporal totalmente aleatoria
print(aleatorio)## Time Series:
## Start = 1
## End = 100
## Frequency = 1
## [1] 0.98774045 2.30305902 2.63641713 0.20642226 1.19894053
## [6] 0.10965624 2.59409187 -0.38686658 0.26106560 1.90999680
## [11] -0.04317301 0.52297030 1.89569828 1.62804745 -1.12651675
## [16] 0.11974763 0.14294881 2.92239597 -1.34521152 0.77074962
## [21] -1.91448657 -0.24371338 0.22218526 0.08032926 0.69886201
## [26] -0.07704617 -0.27697819 -0.13605829 0.01570991 0.72711151
## [31] 0.27850205 0.40099156 -1.66879047 -0.21830735 -1.66006470
## [36] -1.49610220 0.79343682 -1.68216448 0.48056236 -0.18946779
## [41] 0.37695342 0.69819127 -0.15755725 0.80179430 1.62321711
## [46] 1.61495222 -0.26766444 0.30133654 -0.85143769 0.57862783
## [51] -1.29740036 -1.02612585 2.28446087 1.10875990 1.96438933
## [56] -1.03244542 -2.18020840 0.95117935 1.01850255 -1.19561725
## [61] -0.08166370 0.68383541 -0.59202470 -1.73426363 -0.04481183
## [66] 0.53268878 0.39607558 0.36874644 -1.08092582 -0.71681877
## [71] -0.15064003 -1.31195505 0.73892051 0.57235440 -0.87086375
## [76] 1.26895119 -0.40109010 -0.17781604 1.11782754 -1.76333283
## [81] 0.47072571 -0.29399474 -0.01077399 1.21481598 0.46415216
## [86] 0.21240916 0.03113017 -1.14232833 1.02743417 -0.66230460
## [91] 0.33361987 -0.64643847 -0.69601923 1.10313149 0.15473126
## [96] 1.81016203 -0.73860417 -2.33047475 3.03250692 -0.07569329(acf(aleatorio,type='covariance'))
##
## Autocovariances of series 'aleatorio', by lag
##
## 0 1 2 3 4 5 6 7
## 1.29145 -0.06761 0.02910 0.06892 -0.03803 0.25783 -0.03690 0.05633
## 8 9 10 11 12 13 14 15
## 0.16371 -0.06727 0.00939 0.03917 0.02704 0.13165 -0.04394 0.04797
## 16 17 18 19 20
## 0.01884 -0.05720 -0.00505 -0.20113 -0.06752Coeficientes de autocorrelación
También usaremos la función acf() y asumimos una estacionariedad débil. Estimamos los coeficientes de correlación para un intervalo k, como el cociente entre el coeficiente de autocovarianza en el intervalo k y el coeficiente de autocovarianza para k=0, rk = cK/c0. Este coeficiente de autocorrelación siempre está entre -1 y 1. Al grafico de los coeficientes de autocorrelación se le llama correlograma. Siempre empieza con 1, ya que r0=c0/c0=1. Vamos a practicar con los coeficientes de autocorrelación con el proceso aleatori.
(acf(aleatorio,main="Correlograma de un proceso aleatorio")) #correlograma de un proceso aleatorio
##
## Autocorrelations of series 'aleatorio', by lag
##
## 0 1 2 3 4 5 6 7 8 9
## 1.000 -0.052 0.023 0.053 -0.029 0.200 -0.029 0.044 0.127 -0.052
## 10 11 12 13 14 15 16 17 18 19
## 0.007 0.030 0.021 0.102 -0.034 0.037 0.015 -0.044 -0.004 -0.156
## 20
## -0.052Como vemos, r0=1, y luego todos los valores son poco significativos (están por debajo de las lineas azules discontinuas que marcan el nivel de significación). Esto sucede porque nuestro proceso es completamente aleatorio.
Como modelar series temporales. Random Walks y una introducción a medias moviles (MA)
Random Walks
El modelo de Random Walks asume que el valor de una serie temporal en el momento t es el valor de la serie en el momento t-1 más un valor aleatorio (ruido blanco). Básicamente el modelo estima que el valor de una serie en el tiempo t es la suma del ruido blando desde el tiempo 0 hasta el tiempo t. Con este patron, la media va a depender de t, asi como la varianza, y por tanto, este modelo no genera una serie estacionaria. Vamos a simularlo en R.
x <- NULL
x[1] <- 0
for(i in 2:1000){
x[i] <- x[i-1] + rnorm(1)
}
print(x)## [1] 0.0000000000 0.1276700479 0.3978253412 0.5324508035
## [5] 1.4860825253 3.9030266168 5.5304365793 6.2218541811
## [9] 8.0916183709 10.2685315697 10.1808465550 9.2837037500
## [13] 10.8473034780 11.0207893009 11.2489997173 9.8909435862
## [17] 9.1855230716 8.7244469286 8.1178887809 7.4900354067
## [21] 7.6603548521 6.7304707465 6.3020256778 4.8974096804
## [25] 5.6528268482 5.7791218852 5.9868505212 4.8832616223
## [29] 4.5516665777 4.7189309874 3.2165596231 2.7654193952
## [33] 2.6561432274 0.9993517685 0.3618611101 -1.0783502144
## [37] -1.3911526934 -0.6525125207 -0.0365004532 -0.8032103891
## [41] -1.5888980892 -1.1567360776 0.1142677953 1.2607989377
## [45] -0.4442073362 1.4797703336 0.7065726813 0.9637651298
## [49] 2.0762067524 2.4969308884 0.7824731094 1.0095489798
## [53] 1.3184199820 1.4210585414 1.2839940398 0.8292691419
## [57] 3.0194324758 1.8970354589 1.4973524010 2.8973022084
## [61] 2.8653192638 3.9898297347 6.2480981852 6.0854425841
## [65] 3.7941705688 4.7464927007 4.7161302746 3.5996421908
## [69] 4.2718051032 5.0743329085 4.3414748034 4.9426950206
## [73] 4.0963078686 3.2185655606 3.1124707015 2.0497739287
## [77] 1.9187113153 0.0475641032 0.2286733186 0.4624147087
## [81] -1.2613533522 -2.2436983128 -1.8015854327 -1.0355870657
## [85] -0.2435867289 0.0006210528 0.3450631445 1.2641148478
## [89] -0.0068480138 0.4657828660 1.4378865864 2.8338253016
## [93] 1.0553355577 0.5895848847 -0.4623956878 -1.3409152491
## [97] -0.9934468980 -0.6461803312 -0.0840996816 1.2680498255
## [101] 0.8826674341 0.6742739711 2.0893115499 0.7923413206
## [105] 1.8137566632 1.2693088649 2.7939521314 4.1450668587
## [109] 6.3307641165 5.1005652740 4.5657312461 5.5345328288
## [113] 6.5977206717 8.4326211327 9.4901943760 9.3474658531
## [117] 10.5280912699 9.9706603910 9.1274791295 8.4727495084
## [121] 8.4068118104 8.6626588612 8.2594627179 8.2962026098
## [125] 7.4827271921 8.7669514801 8.9058099293 8.8976341233
## [129] 9.4840872236 9.0245564385 10.2840224099 11.0857664105
## [133] 10.1059651443 9.2403131743 8.7217391603 8.6667299153
## [137] 7.8096435255 7.5639468088 7.8583134480 8.1991250396
## [141] 7.9128722216 7.8494048816 8.7495924283 9.1031570838
## [145] 9.8663995840 10.4845344438 11.7525956530 12.1272603990
## [149] 11.9523740862 12.7494051335 13.0361754168 13.8985746596
## [153] 14.2293111493 15.7961617586 16.5928160782 15.0843444670
## [157] 14.2566523280 12.6145082063 11.8569516679 13.1046814617
## [161] 13.7097793819 13.9288871218 13.1493396932 15.0449668533
## [165] 15.5516082130 16.0292022144 17.2077528181 17.4225803932
## [169] 18.3262812341 16.0707605593 16.3522610168 16.4726861366
## [173] 16.7236992553 15.9989861670 15.2010283018 15.7325974001
## [177] 14.1902794164 14.2506889142 13.0684450224 14.5760116873
## [181] 12.7071083878 11.1067003142 11.3191180217 11.1966132947
## [185] 10.4272919474 9.7447175404 10.2580050289 11.8927757111
## [189] 11.3563575679 10.5077381927 11.7294724371 12.1331899185
## [193] 12.6777833780 13.0673888764 12.6452205833 12.8786911706
## [197] 11.8917778163 12.8776046932 12.3682726379 12.4884327705
## [201] 11.5660456858 10.9026060634 11.1242495047 11.9870956366
## [205] 13.0188678647 13.5013690319 13.3047235950 13.0759528714
## [209] 12.9694271977 13.6999788612 14.7242541017 14.2099925014
## [213] 13.2418327114 13.5341098941 12.8489873953 13.2814087715
## [217] 13.9845992738 11.5575238201 11.9331561492 13.0089177828
## [221] 12.3010166775 12.5371488970 13.2328741754 12.6512470143
## [225] 11.6004540782 10.2904939201 9.4366770621 10.0357232979
## [229] 9.7013881478 9.1879880048 9.3171542150 9.0014034571
## [233] 8.9509989411 8.5143843622 9.2611080059 8.1394933389
## [237] 6.8981054174 7.5648455497 6.9996853376 7.4536396789
## [241] 7.0078276370 7.3570395704 8.2942875875 8.0266329674
## [245] 9.7250309335 9.7927710991 11.0616733251 11.6974419843
## [249] 11.7873677333 12.0348803630 11.8664483517 12.0051623426
## [253] 11.7322581386 10.2472323500 10.6013403325 11.5659016424
## [257] 10.6781344965 8.9908032937 9.1373194983 8.5617575458
## [261] 9.0450879374 8.2077725679 9.8907631694 9.8247984276
## [265] 9.7573626121 9.6614768330 10.2308794562 11.0508714661
## [269] 10.4258667289 10.5619313408 11.2246296597 11.5665121551
## [273] 12.4906863469 10.8296662000 10.6700695887 11.6475090717
## [277] 12.8047547489 12.4693567887 11.6926025017 11.1100577988
## [281] 11.4588262134 12.2812220559 13.1125833322 14.1372982408
## [285] 14.3768968851 13.3176551706 13.6364232704 13.2406008584
## [289] 14.1683676396 15.5253404320 16.0995160631 15.7446484523
## [293] 14.7809228620 15.3332587181 14.7924234605 13.7331441479
## [297] 14.6336991177 14.1948802676 14.3348982992 14.4935734063
## [301] 15.4009436559 17.0210697106 17.0677585194 17.7823450432
## [305] 18.8244748710 19.0064976991 19.7344958849 19.7848854045
## [309] 18.3049984913 19.5193675946 18.8804416126 19.0420148751
## [313] 19.5733288491 20.0592700821 20.2637719185 19.8071351218
## [317] 19.1258268620 19.3990204529 19.9467163190 19.9180301276
## [321] 18.1562601483 17.4916538019 17.7122507594 16.8731504733
## [325] 18.1112040303 17.9205418476 17.8213522245 18.6229906482
## [329] 17.6557662391 15.8911794871 16.3644456443 15.2904245766
## [333] 15.4350130825 16.6186820466 16.4800359233 17.2957251629
## [337] 18.0797540762 16.6292227580 15.6100510741 14.1974536588
## [341] 13.3206005896 13.5547751622 13.1183284062 11.4717556663
## [345] 8.7663275023 10.1974937789 9.6850238348 8.5438536249
## [349] 7.3576594377 8.3108027252 7.9874926901 8.5584595769
## [353] 8.4894398689 8.9527475617 10.7210510914 9.6643694048
## [357] 9.6068328654 10.1038130036 10.0363851074 11.8386723246
## [361] 11.9693997978 11.9409175427 10.2793778187 10.7041193204
## [365] 10.4348748852 9.9993329612 10.6122088246 9.5327680530
## [369] 9.4664794314 9.2930908685 9.1192129677 6.5453889380
## [373] 7.1146331359 8.8178777321 8.7618720380 8.6002368586
## [377] 6.8952453823 8.7173425467 8.3201204681 8.7401642628
## [381] 9.6303513114 8.9255948161 8.5794230287 9.0600168025
## [385] 10.7873983482 8.8525797804 8.1721578803 8.6446236433
## [389] 8.8871271729 7.3958851737 6.6758545193 5.7113988470
## [393] 6.5421594269 6.2555882047 9.4237423875 10.7583507922
## [397] 10.5998298266 11.3576039035 10.7704245143 11.0432943482
## [401] 10.8714538194 10.1246306861 8.5045374742 7.9789298289
## [405] 8.8591940851 9.3101468083 11.9787724137 10.0869651168
## [409] 9.7563943571 9.8074670630 10.2456115700 9.5380148425
## [413] 8.9562608680 8.8803153593 10.6753638094 11.7342749023
## [417] 11.9763897402 11.7645456034 10.7940750556 11.7682647159
## [421] 12.5781760441 12.7877481313 12.9029350124 13.0394759669
## [425] 11.8313454915 10.9107674670 11.0140954551 9.8175606386
## [429] 10.2707372270 8.9496996039 10.2785672041 9.3170445526
## [433] 10.3986296948 8.8935127907 8.2947389057 6.7046171698
## [437] 6.6327497649 6.3600081297 6.2359161898 5.3223149089
## [441] 5.6534778660 6.0106836883 5.7836723571 4.6915794694
## [445] 5.6028864173 5.5873738550 6.7156805125 7.8125882730
## [449] 7.7266914555 8.0126921758 9.1124894958 8.5880134202
## [453] 9.0719853411 9.1741185215 10.1513840907 11.7782439009
## [457] 10.3964272037 12.3673756350 12.4908422187 12.0600126982
## [461] 12.7430246785 13.8202303731 13.5779961647 13.6435067337
## [465] 12.9397473715 12.5661578886 11.9395665844 13.0612824601
## [469] 13.6950922215 11.2219819269 12.9734980500 11.6696391569
## [473] 11.5600989857 10.7312720295 11.1000671259 10.4471841105
## [477] 9.7874916718 11.3251317770 11.0427191769 10.8585577492
## [481] 9.3765853042 8.4685733774 8.3279647482 8.6814964971
## [485] 7.9765902902 9.0475429769 7.2945324065 7.2106229103
## [489] 8.0216441017 8.2052298278 8.0952075069 9.5352060056
## [493] 10.3981283604 10.5344192372 8.9715563310 9.9125249351
## [497] 8.0746840335 6.9292251430 5.1220646146 4.6410476983
## [501] 4.7744355709 5.8555594897 6.3711884095 7.6455100930
## [505] 5.8226323336 4.7251948313 3.6102605336 3.8330215119
## [509] 2.5061305595 2.6629409995 2.5156428225 3.2838320646
## [513] 2.4520842990 2.0407988098 1.7885254684 1.8846392944
## [517] 2.1097520636 1.2995096405 2.1958126281 1.7834520609
## [521] 2.6346124623 2.2980822174 1.0838456509 0.9851751447
## [525] 1.8391374856 2.8923251234 1.7962990800 3.3904219324
## [529] 1.4682361575 1.6395873154 1.4161348830 0.4994543491
## [533] 1.1763947868 1.9116208248 2.1627635305 2.6927731572
## [537] 3.8314363485 3.8619171018 2.7948246680 2.8621590156
## [541] 2.9811231191 1.9023610043 1.9138345029 1.0342096472
## [545] 1.1963114768 1.8649959003 2.1563395312 1.9363973380
## [549] 1.8445046504 0.3225946873 -0.8811459158 -2.4979918823
## [553] -2.6725471180 -4.2651789873 -3.1828685042 -3.4796781735
## [557] -4.8132856533 -3.0644567432 -1.9249663732 0.4455308789
## [561] -0.5295681207 0.6209934326 -0.0847822843 0.4409867711
## [565] 0.8437624311 -0.3128598731 -1.5887793427 -1.2339413517
## [569] -1.4214306849 -0.9598664700 -1.8489934417 -1.9221733800
## [573] -1.7292038536 -0.4521437561 -0.3415204168 1.6516726641
## [577] 1.9575924860 2.2788048300 3.2952828775 4.2478024989
## [581] 3.2837149100 4.4007240917 4.2273569613 3.3562713681
## [585] 3.0939361024 4.0508879391 4.5342739658 5.5065639986
## [589] 5.0597802633 4.9949596768 5.6122735971 5.7945137846
## [593] 7.1248033173 7.2273975991 7.2240357730 7.4040402511
## [597] 7.3081238509 7.3703274364 7.6635005327 8.3405502967
## [601] 9.2879008232 7.4999381229 8.4486674504 8.7062222561
## [605] 7.6275921148 7.1436983692 6.9796997059 8.9543063444
## [609] 9.2308174833 7.7357410767 8.8032894519 10.6988985643
## [613] 10.9433520391 10.3289532981 8.2599834973 6.8609422513
## [617] 7.2100056706 8.6506161613 9.3547254372 10.6004281854
## [621] 11.7055531797 13.7520040367 12.6341642157 11.4058825827
## [625] 13.6435141309 15.2333139915 13.8150474360 14.8928432971
## [629] 13.9872928196 12.8160656933 12.7654962955 14.1008845951
## [633] 14.3845946293 15.9128438870 17.1409617130 16.5540341146
## [637] 16.1535377730 17.4248583903 17.5816332241 17.9048545141
## [641] 18.5911124566 17.7977296992 17.5772707116 16.9245229012
## [645] 17.4937154906 17.9347027601 17.0995466474 16.3696858622
## [649] 15.4754057585 13.4993848334 12.7311215499 13.2839973360
## [653] 13.2750545104 13.3200873052 13.7019244624 12.1202667374
## [657] 9.6992929811 9.0086052347 7.0367737927 7.8626622067
## [661] 6.9266597751 8.1346060726 6.2535892993 6.9319160392
## [665] 5.9254659353 3.0546401652 2.7048413726 2.0102695500
## [669] 0.9019385289 0.6577399871 0.8465676334 1.8214205986
## [673] 1.5147742699 1.0411812435 0.7526006910 -0.1781782853
## [677] 0.0065881206 -0.2811207114 -1.3727313729 -1.1836570550
## [681] -2.0821807632 -1.9750670480 -1.7838309915 -3.0047516050
## [685] -2.3358102155 -2.5162787807 -3.3797513271 -4.1010825615
## [689] -5.3475443512 -4.6784337069 -5.3902730945 -3.9106302833
## [693] -3.7872210759 -3.8073781884 -4.1646136345 -4.1712875056
## [697] -3.9541236795 -5.8922471500 -4.1215298797 -2.3633326540
## [701] -5.2700204588 -7.0242429783 -5.9985431216 -5.8282031226
## [705] -4.8903903959 -4.9587122617 -3.6524048867 -3.3870941812
## [709] -1.5064997093 -3.3620974310 -3.6816573492 -4.9082413335
## [713] -5.4469773236 -6.3919546601 -5.6429917828 -6.0740963434
## [717] -4.0597588104 -2.9755183589 -2.5786841390 -2.3994993358
## [721] -3.4980548977 -3.2173795277 -2.1951199412 -2.4585964409
## [725] -2.8701834407 -4.6954814626 -4.9799435126 -5.3194307299
## [729] -7.1765964021 -6.6772494711 -6.2713967813 -4.2625094298
## [733] -3.9956749355 -3.1177407273 -3.2306382582 -4.9266552372
## [737] -2.4566546840 -3.2709790516 -4.0742774358 -3.2328091876
## [741] -2.3434884311 -3.4343127532 -4.4132355265 -4.2579664484
## [745] -3.5420700375 -2.0911030685 -1.8263875845 -3.7500412712
## [749] -2.8880539836 -1.8508858614 -0.6630156138 0.9171894381
## [753] -0.6895893659 -0.0639963936 -0.0320136469 -1.5460351682
## [757] -2.2793826118 -0.8827377534 -0.3427213303 -1.3330796353
## [761] -1.7159382077 -2.4463941235 -4.0090912495 -4.2004276102
## [765] -4.9782452971 -5.1440135811 -3.7334152426 -2.6160754550
## [769] -3.4491094708 -1.5469463640 -0.3922898678 -0.5230134080
## [773] 1.2993256489 0.8152909502 1.0276340445 1.4293738897
## [777] 0.5494532759 0.9313754416 1.1159632180 1.4607242012
## [781] 2.1462560244 2.4313923311 3.0879287589 1.5406589380
## [785] -0.7591124586 -0.7316412024 -1.2027965353 -0.5349840572
## [789] 0.0362788646 0.1665376311 -1.4935641846 -0.8593344598
## [793] -0.6451879336 0.5083802694 1.3004352464 1.4954584245
## [797] 3.5657079155 1.7563021208 3.3412937591 2.1746887984
## [801] 0.4194323378 0.3331725971 2.0238987568 0.7290114083
## [805] -0.4199831906 -1.1506563130 -1.7479903154 0.4573235680
## [809] -0.3257500329 0.1648949085 1.0015722748 1.7635339737
## [813] 3.0166911087 2.7671716199 3.8587745499 6.0793630652
## [817] 8.0159049920 8.1267759293 7.2725438541 6.4971823271
## [821] 4.9727135310 3.5330474251 2.3855541548 3.2050944841
## [825] 3.2748290892 3.7436847380 4.1020969724 3.6396105073
## [829] 4.7964975797 5.7577655967 6.4970617648 4.5417829750
## [833] 4.2090651405 5.1464378843 6.3797994838 5.5184272527
## [837] 5.2483692605 5.4420357784 5.5451741900 6.9786206819
## [841] 7.1631912615 8.3894660627 8.3084244586 8.2777139054
## [845] 9.1587074151 9.5037432629 7.9226729565 7.0394707310
## [849] 6.9090248040 6.8200761622 7.0614997194 7.1272158281
## [853] 6.9433332462 7.1636455467 7.6930796035 6.8411535669
## [857] 7.1049622122 6.4929086630 6.3453214638 5.3745282214
## [861] 6.9451310516 7.7274889372 7.5373394181 5.7235527415
## [865] 6.7932615479 8.2774299770 7.9040086127 8.6663288215
## [869] 9.9065278919 10.0388643891 9.4676452031 9.9730492799
## [873] 9.2512309760 10.4372068641 10.0338116522 10.2075919632
## [877] 12.5159475754 12.6563058685 15.4894559612 14.9224127612
## [881] 14.5951352212 14.6535746314 14.6938426390 14.4001311840
## [885] 14.2419528017 13.5773723220 14.6328199215 13.6706403367
## [889] 13.6246077094 13.1616752298 14.2404780840 14.4177543867
## [893] 14.3968489411 15.1933584263 14.4588921519 15.1383594219
## [897] 15.8493438587 17.4312991031 16.0717800992 14.3235756680
## [901] 15.4384257877 14.2863031756 13.4638900081 13.3261784052
## [905] 12.2135024586 11.3898673269 11.8356101953 11.3878470377
## [909] 11.0246327253 11.4402752629 11.2373545406 10.9399688303
## [913] 12.0262002702 12.0254331404 12.3852038688 12.2327269327
## [917] 13.1688924269 11.9025242829 12.0338943849 11.2963498474
## [921] 11.8838720034 10.9866367527 11.0957154383 12.2292228330
## [925] 12.8971037470 12.0919537801 12.2556223612 13.3016267893
## [929] 13.6420880009 13.6854671051 11.0618116552 13.7645377982
## [933] 14.7930214870 14.3258978698 12.9265090712 13.1790340691
## [937] 13.4447606973 11.3246787886 9.9385960540 10.4245978829
## [941] 8.0125116487 7.4641302046 7.7728354192 8.6657983340
## [945] 9.4206666580 9.5900131709 9.6664592243 8.8896206989
## [949] 9.7529976914 10.5685764534 9.4194086351 9.0989673937
## [953] 7.6878765465 6.7330390822 7.2506128578 7.1865818085
## [957] 8.6960945990 7.6042586034 7.3921674859 6.2964577237
## [961] 5.3081490415 5.9999563046 6.7645677164 5.1492026426
## [965] 6.0547679937 5.6887148382 6.7846607459 8.8155825221
## [969] 9.2820080835 9.9432159025 8.9001391314 8.6274939359
## [973] 8.4537857313 9.5203523242 8.6531433591 8.4014411889
## [977] 8.4545204666 9.4414897210 9.6904963558 10.5164588128
## [981] 10.3022811740 11.2579950611 11.2348289081 10.8739006155
## [985] 11.9606778983 11.3310975309 9.8811718915 10.7695891326
## [989] 10.6170268932 13.0361852624 11.2391200641 9.4250908007
## [993] 9.0084164685 8.0717081554 8.6794078013 8.0656373447
## [997] 9.4767380296 11.6819600407 11.5846367141 12.5045628161#le damos formato de serie temporal
random_walk <- ts(x)
#dibujamos la serie temporal, que se verá que no es estacionaria
plot(random_walk,main="Random Walk",xlab="Dias",ylab=" ",col="blue",lwd=2)
#aunque los coeficientes de autocorrelación se definen para series estacionarias, vamos a probarlo con el random walk
acf(random_walk)
Normalmente queremos trabajar con series temporales estacionarias, al menos de forma débil, y en la realidad estas no suelen existir. Vamos a ver como transformamos una serie en estacionaria, eliminando la tendencia. Si hacemos la diferencia de la serie temporal consigo misma, pero desplazandola una unidad de t, nos quedará que Xt - X(t-1) = Z(t), y como sabemos Z es un ruido blanco, y por tanto, estacionario. Vamos a verlo en R.
plot(diff(random_walk)) #la función diff sin parametros nos da las diferencias entre valores t contiguos
#calculamos los coeficientes de autocorrelación de la diferencia
acf(diff(random_walk))
Introducción a las medias móviles (MA)
Un modelo de medias moviles de orden q (MA(q)) es un modelo en el cual el valor actual de la variable X(t) es influeciado por los valores de q valores anteriores de X (o del ruido blanco que la genera, si tenemos una serie de estilo random walk), cada uno con un coeficiente distinto.
# Generate noise
noise=rnorm(10000)
# Introduce a variable
ma_2=NULL
# Loop for generating MA(2) process
for(i in 3:10000){
ma_2[i]=noise[i]+0.7*noise[i-1]+0.2*noise[i-2]
}
# Shift data to left by 2 units
moving_average_process=ma_2[3:10000]
# Put time series structure on a vanilla data
moving_average_process=ts(moving_average_process)
# Partition output graphics as a multi frame of 2 rows and 1 column
par(mfrow=c(2,1))
# plot the process and plot its ACF
plot(moving_average_process, main='A moving average process of order 2', ylab=' ', col='blue')
acf(moving_average_process, main='Correlogram of a moving average process of order 2')
Como vemos, un proceso MA(2) muestra en su ACF valores significativos hasta la diferencia 2. A partir de ahí ya no hay valores significativos. En general, un proceso MA(q) tiene valores significativos en su ACF hasta el valor q de diferencia.
Semana 3 - Estacionariedad, procesos MA(q) y AR (p)
Estacionariedad: generalizando de lo individual al grupo
Estacionariedad: intuición y definición
Un proceso estocástico es una familia de variables aleatorias estructuradas con una base temporal, y conocemos la relación entre las variables en distintos puntos temporales. Pueden ser discretos o continuos. Los procesos estocáticos con realmente complejos de manejar.
Normalmente tenemos una secuencia de datos observados (realizaciones) e intentamos averiguar las propiedades del proceso que los ha generado de esta simple trayectoria.
¿Como podemos inferir las propiedades del proceso estocástico de una simple realización?
Un proceso estrictamente estacionario si la distribución de un conjunto de realizaciones es igual de a otro conjunto en un tiempo distinto. Es decir, las propiedades estadisticas del proceso estocástico son independientes del tiempo. Que estén igualmente distribuidas no significa que sean independientes. Además, la media y la varianza del proceso estocástico será constante. También la función de autocovarianza solo dependerá de la distancia en el tiempo de ambas variables, no del momento temporal concreto.
Se define la estacionariedad débil cuando solo la media es constante y si la autocovarianza solo depende de la separación temporal de las variables que forman el proceso estocástico. Esto implica también una varianza constante.
Estacionariedad: primeros ejemplos - ruido blanco, random walks
El ruido blanco es estacionario, de media cero y varianza constante. Por contra, los random walks no son estacionarios, ya que su media será la media de las variables aleatorias que añadimos por el t, asi como la varianza será la varianza de la variable aleatoria por t. Por tanto, tanto media como varianza dependen del tiempo t. Los procesos de medias moviles MA(q) son estacionarios.
Estacionariedad: primeros ejemplos - ACF de las medias moviles MA(q)
Recordemos que un proceso de medias moviles de orden q, MA(q), se define como X(t))=B0Z(t)+B1Z(t-1)+…+BqZ(t-q)*, donde Z(t) es una variable independientemente distribuida de media 0 y varianza sigma2. El orden q nos dice lo lejos que llega la influencia de la variable aleatoria Z en cada momento. Este proceso estocástico es estacionario, ya que la media no depende del tiempo, y la autocovarianza solo depende de la separación temporal de las variables aleatorias.
Backward shift operator para procesos MA(q) y AR(p)
Series y su representación
Secuencia - lista de numeros. Si existe el limite de la secuencia, se dice que la secuencia converge.
Sumas parciales - suma de algunos elementos consecutivos de una secuencia. Estas sumas parciales forman una nueva secuencia llamada serie. Si esta secuencia de sumas parciales converge a un numero, diremos que es convergente, sino que es divergente.
Una serie es absolutamente convergente si la suma de los valores absolutos de los miembros de la serie es convergente.
Un caso importante de series son las series geométricas. Si el parametro r es menor de 1, la serie geométrica será convergente. Ciertas funciones se pueden representar como series geométricas, siempre que el valor absoluto de r sea menor de 1.
Backward shift operator
Si tenemos un proceso estocástico X(t) - X1, X2,X3,X4…
El operador Backward shift se define como BX(t) = X(t-1). En general, Bk X(t) = X(t-k).
Introducción a la invertibilidad
Un proceso MA(1) se puede convertir, o invertir hacia atras, a un proceso AR(inf).
Un proceso invertible asegura que el ACF observado corresponde a un solo proceso MA(q).
Dualidad
La dualidad se refiere a la relación entre los procesos MA y AR.
Si un proceso MA(q) es invertible, entonces se puede expresar como un modelo AR(inf).
Si un proceso AR(p) es estacionario, entonces se puede expresar como un modelo MA(inf).
Procesos autoregresivos AR(p)
Definición, simulación y primeros ejemplos
Un proceso autorregresivo AR(p) es aquel en el que el valor actual de la serie se forma con la suma ponderada de p valores anteriores en el tiempo más un termino de ruido blanco Zt. Los procesos autorregresivos no han de ser necesariamente estacionarios. Vamos a simular un proceso AR(1) y veremos que dependiendo del valor del coeficiente el proceso puede ser o no estacionario.
set.seed(2016)
N=1000
#fijamos el coeficiente a 0.4
phi=0.4
Z <- rnorm(N,0,1)
X <- NULL
X[1] <- Z[1]
for(t in 2:N){
X[t] <- Z[t] + phi*X[t-1]
}
X.ts <- ts(X)
par(mfrow=c(2,1))
plot(X.ts,main="AR(1) serie temporal con ruido blando, phi=0.4")
acf(X.ts,main="AR(1) serie temporal con ruido blando, phi=0.4")
par(mfrow=c(1,1))
#fijamos ahora el coeficiente a 1
phi=1
set.seed(2016)
Z <- rnorm(N,0,1)
X <- NULL
X[1] <- Z[1]
for(t in 2:N){
X[t] <- Z[t] + phi*X[t-1]
}
X.ts <- ts(X)
par(mfrow=c(2,1))
plot(X.ts,main="AR(1) serie temporal con ruido blando, phi=1")
acf(X.ts,main="AR(1) serie temporal con ruido blando, phi=1")
par(mfrow=c(1,1))Ahora vamos a simular un proceso AR(2).
set.seed(2017)
X.ts <- arima.sim(list(ar=c(0.7,0.2)),n=1000)
par(mfrow=c(2,1))
plot(X.ts,main="AR(2) serie temporal con phi1=0.7, phi2=0.2")
acf(X.ts, main="Autocorrelación de AR(2)")
Vemos que en el proceso AR(2) la autocorrelación tarda tiempo en desaparecer, que hace que el proceso no sea del todo estacionario. Para que un proceso AR(2) sea estacionario, se debe cumplir:
- -1 < phi2 < 1
- phi2 < 1 + phi1
- phi2 < 1 - phi1
Vamos a cambiar los coeficientes del proceso AR(2) a 0.5 y -0.4. Vemos que la autocorrelación baja mucho más rapidamente que en el caso anterior.
set.seed(2017)
phi1=0.5
phi2=-0.4
X.ts <- arima.sim(list(ar=c(phi1,phi2)),n=1000)
par(mfrow=c(2,1))
plot(X.ts,main=paste("AR(2) serie temporal con phi1=",phi1,"-phi2=",phi2))
acf(X.ts, main="Autocorrelación de AR(2)")
Procesos AR, backward shift operator y ACF
En los procesos AR, la media es cero, y la varianza es la varianza del ruido blanco multiplicada por la suma infinita de los coeficientes del proceso MA equivalente. PAra que el proceso AR sea estacionario, esta suma infinita (es decir, la varianza) debe converger.
Introduccion a las ecuaciones Yule-Walker
Ecuaciones diferencia
Una ecuación diferencia es aquella que contiene algún termino recursivo. Por ejemplo, an=5a(n-1)-6a(n-2). Estas ecuaciones se resuelven utilizando un cambio de varible: an=lambda^n. Con esto nos permite tener una expresión analítica para una serie cualquiera.
Ecuaciones Yule-Walker
Utilizaremos las ecuaciones de Yule-Walker para obtener el ACF de un proceso autorregresivo AR.
Semana 4 - procesos AR(p), ecuaciones Yule-Walker y PACF
Employ PACF to estimate the order of AR( p ) processes.
Autocorrelacion parcial (PACF) - primeros ejemplos
Para los procesos de medias moviles de orden q sabiamos que si el ACF tiene 3 valores significativos, el proceso sera de orden 3, MA(3). Podemos hacer lo mismo para procesos autorregresivos? PAra ello tenemos el PACF (Partial ACF). Veamos algún ejemplo donde ya se puede intuir algo:
par(mfrow=c(3,1))
phi.1 = .6
phi.2 = .2
data.ts = arima.sim(n = 500, list(ar = c(phi.1, phi.2)))
plot(data.ts, main=paste("Autoregressive Process with phi1=",phi.1," phi2=",phi.2))
acf(data.ts, main="Autocorrelation Function")
acf(data.ts, type="partial", main="Partial Autocorrelation Function") #vemos 2 valores significativos
##ahora un proceso AR de orden 3
par(mfrow=c(3,1))
phi.1 = .9
phi.2 = -.6
phi.3 = .3
data.ts = arima.sim(n = 500, list(ar = c(phi.1, phi.2, phi.3)))
plot(data.ts, main= paste("Autoregressive Process with phi1=",phi.1, "phi2=",phi.2," phi3=",phi.3))
acf(data.ts, main="Autocorrelation Function")
acf(data.ts, type="partial", main="Partial Autocorrelation Function") #vemos 3 valores significativos
Ahora vamos a trabajar con la serie temporal de los precios del trigo para cerveza de “Weather and Harvest Cycles” (Beveridge, 1921). Cargamos las serie temporal y hacemos algunos gráficos:
beveridge = read.table("beveridge.txt", header=TRUE)
beveridge.ts = ts(beveridge[,2], start=1500)
#evolución de los precios del trigo
plot( beveridge.ts, ylab="price", main="Beveridge Wheat Price Data")
beveridge.MA = stats::filter(beveridge.ts, rep(1/31, 31), sides = 2)
lines(beveridge.MA, col="red")
#Vamos a escalar los valores de la serie temporal por su media
par(mfrow=c(3,1))
Y = beveridge.ts/beveridge.MA #con este escalado intentamos hacer el proceso más estacionario
plot( Y, ylab="scaled price", main="Transformed Beveridge Wheat Price Data")
acf(na.omit(Y),main="Autocorrelation Function of Transformed Beveridge Data")
acf(na.omit(Y), type="partial",main="Partial Autocorrelation Function of Transformed Beveridge Data") #
#Podemos usar la función ar() para estimar los coeficientes del proceso AR de la serie escalada
ar(na.omit(Y), order.max = 5) #estima solo dos coeficientes##
## Call:
## ar(x = na.omit(Y), order.max = 5)
##
## Coefficients:
## 1 2
## 0.7232 -0.2949
##
## Order selected 2 sigma^2 estimated as 0.027Vemos que la funcion del proceso AR solo devuelve 2 coeficientes. y vemos que el PACF solo tiene dos valores significativos. Por lo tanto, podemos decir que el orden p de un proceso autorregresivo AR(p) es igual al numero de valores significativos del PACF.
Autocorrelacion parcial y concepto PACF
Vamos a usar el dataset bodyfat de la libreria isdals. Este set de datos tiene 4 variables: la grasa corporal, la anchura del triceps, la circunferencia del muslo y la circunferencia del antebrazo. Vamos a ver, sin que sea una sorpresa, que las variables están relacionadas.
data(bodyfat)
attach(bodyfat)
pairs( cbind( Fat, Triceps, Thigh, Midarm) ) 
cor(cbind( Fat, Triceps, Thigh, Midarm))#vemos que Fat, Triceps y Thigh están altamente correladas## Fat Triceps Thigh Midarm
## Fat 1.0000000 0.8432654 0.8780896 0.1424440
## Triceps 0.8432654 1.0000000 0.9238425 0.4577772
## Thigh 0.8780896 0.9238425 1.0000000 0.0846675
## Midarm 0.1424440 0.4577772 0.0846675 1.0000000#vamos a modelar las variables Fat y Triceps en función de Thigh, al estar muy relacionadas es posible con una simple regresión lineal
Fat.hat = predict(lm(Fat~Thigh))
Triceps.hat = predict( lm(Triceps~Thigh) )
cor( (Fat- Fat.hat), (Triceps- Triceps.hat) ) #returns 0.1749822## [1] 0.1749822#podemos controlar Far y Triceps sacando Thigh
#Que pasaría si ahora la estimación la hacemos con Thigh y Midarm
Fat.hat = predict(lm(Fat~Thigh+Midarm))
Triceps.hat = predict( lm(Triceps~Thigh+Midarm) )
cor( (Fat- Fat.hat), (Triceps- Triceps.hat) ) #returns 0.33815## [1] 0.33815pcor( cbind( Fat, Triceps, Thigh, Midarm) ) #confirma la correlación## $estimate
## Fat Triceps Thigh Midarm
## Fat 1.0000000 0.3381500 -0.2665991 -0.3240520
## Triceps 0.3381500 1.0000000 0.9963725 0.9955918
## Thigh -0.2665991 0.9963725 1.0000000 -0.9926612
## Midarm -0.3240520 0.9955918 -0.9926612 1.0000000
##
## $p.value
## Fat Triceps Thigh
## Fat 0.0000000 0.169911065182822285102304 0.284894370161702492616485
## Triceps 0.1699111 0.000000000000000000000000 0.000000000000000001490492
## Thigh 0.2848944 0.000000000000000001490492 0.000000000000000000000000
## Midarm 0.1895628 0.000000000000000007071386 0.000000000000000413417781
## Midarm
## Fat 0.189562846639528437275857
## Triceps 0.000000000000000007071386
## Thigh 0.000000000000000413417781
## Midarm 0.000000000000000000000000
##
## $statistic
## Fat Triceps Thigh Midarm
## Fat 0.000000 1.437266 -1.106441 -1.370142
## Triceps 1.437266 0.000000 46.833317 42.459264
## Thigh -1.106441 46.833317 0.000000 -32.834695
## Midarm -1.370142 42.459264 -32.834695 0.000000
##
## $n
## [1] 20
##
## $gp
## [1] 2
##
## $method
## [1] "pearson"PAra series temporales que se pueden modelar con procesos AR(p) es similar ya que en el dataset anterior teniamos también una relación entre variables datas por un desarrollo polinómico con ciertos coeficientes.
La idea básica de la autocorrelación parcial es encontrar la correlación entre las variables aleatorias Xt y X(t-k), despues de eliminar el efecto de combinación lineal entre las variables aleatorias. Al hacer la resta de la variable (o serie temporal) original con su estimación, eliminamos este efecto lineal. La correlación entre las diferencias de las variables originales con su estimación (al eliminar el efecto lineal) es la correlación parcial (PACF).
Vamos a producir ahora un PACF con una serie temporal modelada como un proceso AR y con arima.sim. Al ser un proceso AR(2), el PACF tiene claramente dos valores significativos.
phi.1 = .6;
phi.2 = -.6;
data.ts = arima.sim(n = 1000, list(ar = c(phi.1, phi.2)))
acf(data.ts, type="partial", main=paste("PACF of Time Series Data, phi1=",phi.1,", phi.2=",phi.2) )
Ecuaciones de Yule-Walker en forma matricial y estimacion de parametros de un proceso AR(p)
Ecuaciones de Yule-Walker en forma matricial
Recordemos que un proceso AR(p) es una proceso que calcula el valor actual en base a una combinación lineal de p valores anteriores más un termino aleatorio de ruido blanco.
También asumimos que el proceso AR(p) es estacionario.
Podemos decir, derivada de las ecuaciones de Yule-Walker, que b=Rphi, donde R es la matriz de estimaciones de la correlación que podemos calcular con ACF, b es el vector de autocorrelación y phi es el vector de coeficientes del modelo AR. La dimensión de R es de pxp. Entonces, para calcular los coeficientes, solo debemos calcular phi=R^(-1)b.
Yule Walker Estimation - AR(2) Simulation
Proceso AR(2) –> xt=phi1x(t-1)+phi2x_(t-2)+z_t z_t~ N(0, sigma^2)
set.seed(2017)
#fijamos los coeficientes
sigma <- 4
phi <- NULL
phi[1:2] <- c(1/3,1/2)
phi## [1] 0.3333333 0.5000000n <- 10000 #numero de puntos
#simulamos el proceso con ARIMA
ar.process <- arima.sim(n,model=list(ar=c(1/3,1/2)), sd=4)
ar.process[1:5]## [1] 4.087685 5.598492 3.019295 2.442354 5.398302#calculamos las estimacinoes de la autocorrelacion
r <- NULL
r[1:2] <- acf(ar.process, plot=F)$acf[2:3]
r## [1] 0.6814103 0.7255825#creamos la matriz R
R <- matrix(1,2,2) # matrix of dimension 2 by 2, with entries all 1's.
R## [,1] [,2]
## [1,] 1 1
## [2,] 1 1R[1,2] <- r[1] # only main diagonal entries are kept
R[2,1] <- r[1] # only main diagonal entries are kept
R## [,1] [,2]
## [1,] 1.0000000 0.6814103
## [2,] 0.6814103 1.0000000#ahora creamos el vector b
b <- matrix(r,nrow=2,ncol=1)# b- column vector with no entries
b## [,1]
## [1,] 0.6814103
## [2,] 0.7255825#resolvemos las ecuaciones de Yule-Walker
#phi.hat <- matrix(c(solve(R,b)[1,1], solve(R,b)[2,1]),2,1)
phi.hat <- solve(R,b)
phi.hat## [,1]
## [1,] 0.3490720
## [2,] 0.4877212#estimamos la varianza
c0=acf(ar.process, type='covariance', plot=F)$acf[1]
var.hat=c0*(1-sum(phi.hat*r))
var.hat## [1] 16.37169#hacemos un grafico con la serie, acf y pacf
par(mfrow=c(3,1))
plot(ar.process, main='Simulated AR(2)')
acf(ar.process, main='ACF')
pacf(ar.process, main='PACF')
Yule Walker Estimation - AR(3) Simulation
Proceso AR(3)
x_t=phi1x_(t-1)+phi2 x_(t-2)+_3*x_(t-3)+z_t z_t~ N(0, sigma^2)
#establecemos los parámetros
set.seed(2017)
sigma <- 4
phi <- NULL
phi[1:3] <- c(1/3,1/2,7/100)
n <- 100000
#simulamos el proceso, generamos la matriz R y el vector b
ar3.process <- arima.sim(n,model=list(ar=c(1/3,1/2, 7/100)), sd=4)
r <- NULL
r[1:3] <- acf(ar3.process, plot=F)$acf[2:4]
r## [1] 0.7859646 0.8180901 0.7369167R <- matrix(1,3,3)
R[1,2] <- r[1]
R[1,3] <- r[2]
R[2,1] <- r[1]
R[2,3] <- r[1]
R[3,1] <- r[2]
R[3,2] <- r[1]
R## [,1] [,2] [,3]
## [1,] 1.0000000 0.7859646 0.8180901
## [2,] 0.7859646 1.0000000 0.7859646
## [3,] 0.8180901 0.7859646 1.0000000b <- matrix(1,3,1)
b[1,1]=r[1]
b[2,1]=r[2]
b[3,1]=r[3]
b## [,1]
## [1,] 0.7859646
## [2,] 0.8180901
## [3,] 0.7369167#resolvemos las ecuaciones de yule-walker
phi.hat <- solve(R,b)
phi.hat## [,1]
## [1,] 0.33812448
## [2,] 0.49849991
## [3,] 0.06849712# sigme estimation
c0 <- acf(ar3.process, type='covariance', plot=F)$acf[1]
var.hat=c0*(1-sum(phi.hat*r))
var.hat## [1] 15.979#hacemos las graficas del proceso
par(mfrow=c(3,1))
plot(ar3.process, main='Simulated AR(3)')
acf(ar3.process, main='ACF')
pacf(ar3.process, main='PACF')
Procesos AR(p) - Ejemplos orientados a datos
Modeling recruitment time series - ajuste de modelo
Vamos a practicar con el set de datos de recruitment que tenemos en el paquete astsa.
Usaremos el principio de parsimonia durante este tipo de analisis:
- Elegir la explicación más simple que encaja con las evidencias.
- La teoria más simple de las consideradas será la preferida.
En el siguiente caso, veremos que según el PACF tenemos un proceso AR(2), e intentariemos ajustar a un proceso de media cero para poder usar las ecuaciones de Yule-Walker en forma de matriz.
my.data <- rec
# Plot rec, acf and pacf
par(mfrow=c(3,1))
plot(rec, main='Recruitment time series', col='blue', lwd=3)
acf(my.data,main="ACF from recruitment time series")
acf(my.data,type="partial",main="PACF from recruitment time series")
# subtract mean to get a time series with mean zero
ar.process <- my.data-mean(my.data)
# ACF and PACF of the process
par(mfrow=c(2,1))
acf(ar.process, main='Recruitment', col='red', lwd=3)
pacf(ar.process, main='Recruitment', col='green', lwd=3) #2 significative values, posible AR(2)
# order
p <- 2
# sample autocorreleation function r
r <- NULL
r[1:p] <- acf(ar.process, plot=F)$acf[2:(p+1)]
cat('r=',r,'\n')## r= 0.9218042 0.7829182# matrix R
R <- matrix(1,p,p) # matrix of dimension 2 by 2, with entries all 1's.
# define non-diagonal entires of R
for(i in 1:p){
for(j in 1:p){
if(i!=j)
R[i,j] <- r[abs(i-j)]
}
}
cat('Matriz R de estimaciones de coeficientes de correlación.\n')## Matriz R de estimaciones de coeficientes de correlación.R## [,1] [,2]
## [1,] 1.0000000 0.9218042
## [2,] 0.9218042 1.0000000# b-column vector on the right
b <- NULL
b <- matrix(r,p,1)# b- column vector with no entries
cat('Vector b de estimaciones de coeficientes de correlación.\n')## Vector b de estimaciones de coeficientes de correlación.b## [,1]
## [1,] 0.9218042
## [2,] 0.7829182# solve(R,b) solves Rx=b, and gives x=R^(-1)b vector
phi.hat <- NULL
phi.hat <- solve(R,b)[,1]
cat('Estimación de los coeficientes del modelo AR.\n')## Estimación de los coeficientes del modelo AR.phi.hat## [1] 1.3315874 -0.4445447#variance estimation using Yule-Walker Estimator
c0 <- acf(ar.process, type='covariance', plot=F)$acf[1]
cat('Estimación del coeficiente 0 de la autocovarianza.\n')## Estimación del coeficiente 0 de la autocovarianza.c0## [1] 780.991var.hat <- c0*(1-sum(phi.hat*r))
cat('Valor de la varianza del ruido blando.\n')## Valor de la varianza del ruido blando.var.hat## [1] 94.17131# constant term in the model
phi0.hat <- mean(my.data)*(1-sum(phi.hat))
cat('Estimación de los coeficientes del modelo AR.\n')## Estimación de los coeficientes del modelo AR.phi0.hat## [1] 7.033036cat('Modelo ajustado con las ecuaciones Yule-Walker para un proceso AR(',p,').\n')## Modelo ajustado con las ecuaciones Yule-Walker para un proceso AR( 2 ).cat("Constant:", phi0.hat," Coefficients:", phi.hat, " and Variance:", var.hat, '\n')## Constant: 7.033036 Coefficients: 1.331587 -0.4445447 and Variance: 94.17131Johnson & Johnson-model fitting
VAmos a ajustar un modelo con los datos de J&J. Este será un dataset que no es estacionario ni en media ni en varianza. Para eliminar este problema, hacemos la transformación log.return, que es tomar logaritmos de los valores y hacer la diferencia entre el valor del momento actual y el anterior. De esta manera, después de aplicar la transformación para tener media cero, podemos aplicar las ecuaciones Yule-Walker de forma matricial.
# Time plot for Johnson&Johnson
plot(JohnsonJohnson, main='Johnson&Johnosn earnings per share', col='blue', lwd=3)
# log-return of Johnson&Johnson
jj.log.return=diff(log(JohnsonJohnson))
plot(jj.log.return, main='Johnson&Johnosn earnings per share log difference', col='red', lwd=3)
jj.log.return.mean.zero=jj.log.return-mean(jj.log.return)
plot(jj.log.return.mean.zero, main='Johnson&Johnosn earnings per share log difference mean zero', col='green', lwd=3)
# Plots for log-returns
par(mfrow=c(3,1))
plot(jj.log.return.mean.zero, main='Log-return (mean zero) of Johnson&Johnosn earnings per share')
acf(jj.log.return.mean.zero, main='ACF')
pacf(jj.log.return.mean.zero, main='PACF') #4 valores significativos en el PACF
# Order
p=4
# sample autocorreleation function r
r=NULL
r[1:p]=acf(jj.log.return.mean.zero, plot=F)$acf[2:(p+1)]
r## [1] -0.50681760 0.06710084 -0.40283604 0.73144780# matrix R
R=matrix(1,p,p) # matrix of dimension 4 by 4, with entries all 1's.
# define non-diagonal entires of R
for(i in 1:p){
for(j in 1:p){
if(i!=j)
R[i,j]=r[abs(i-j)]
}
}
cat('Matriz R de estimaciones de coeficientes de correlación.\n')## Matriz R de estimaciones de coeficientes de correlación.R## [,1] [,2] [,3] [,4]
## [1,] 1.00000000 -0.50681760 0.06710084 -0.40283604
## [2,] -0.50681760 1.00000000 -0.50681760 0.06710084
## [3,] 0.06710084 -0.50681760 1.00000000 -0.50681760
## [4,] -0.40283604 0.06710084 -0.50681760 1.00000000# b-column vector on the right
b=matrix(r,p,1)# b- column vector with no entries
cat('Vector b de estimaciones de coeficientes de correlación.\n')## Vector b de estimaciones de coeficientes de correlación.b## [,1]
## [1,] -0.50681760
## [2,] 0.06710084
## [3,] -0.40283604
## [4,] 0.73144780phi.hat=solve(R,b)[,1]
cat('Estimación de los coeficientes del modelo AR.\n')## Estimación de los coeficientes del modelo AR.phi.hat## [1] -0.6293492 -0.5171526 -0.4883374 0.2651266# Variance estimation using Yule-Walker Estimator
c0=acf(jj.log.return.mean.zero, type='covariance', plot=F)$acf[1]
cat('Estimación del coeficiente 0 de la autocovarianza.\n')## Estimación del coeficiente 0 de la autocovarianza.c0## [1] 0.04365692var.hat=c0*(1-sum(phi.hat*r))
cat('Valor de la varianza del ruido blando.\n')## Valor de la varianza del ruido blando.var.hat## [1] 0.01419242# Constant term in the model
phi0.hat=mean(jj.log.return)*(1-sum(phi.hat))
cat('Estimación del coeficiente 0 del modelo AR.\n')## Estimación del coeficiente 0 del modelo AR.phi0.hat## [1] 0.079781cat('Modelo ajustado con las ecuaciones Yule-Walker para un proceso AR(',p,').\n')## Modelo ajustado con las ecuaciones Yule-Walker para un proceso AR( 4 ).cat("Constant:", phi0.hat," Coefficients:", phi.hat, " and Variance:", var.hat, '\n')## Constant: 0.079781 Coefficients: -0.6293492 -0.5171526 -0.4883374 0.2651266 and Variance: 0.01419242Modelo de la serie LakeHuron
VAmos a ajustar un modelo para la serie temporal del LakeHuron, que mide el nivel del lago (en pies) durante 5 años.
head(LakeHuron)## Time Series:
## Start = 1875
## End = 1880
## Frequency = 1
## [1] 580.38 581.86 580.97 580.80 579.79 580.39plot(LakeHuron)
#la serie tiene cierta tendencia
#eliminamos la tendencia con la diferencia con el valor anterior
df <- LakeHuron
ar.process.LH <- diff(df)
plot(ar.process.LH)
# ACF and PACF of the process
par(mfrow=c(2,1))
acf(ar.process.LH, main='Lake Huron ACF', col='red', lwd=3)
pacf(ar.process.LH, main='Lake Huron PACF', col='green', lwd=3) #2 significative values, posible AR(2)
#tres primeros coeficientes de autocorrelación
acf(ar.process.LH,plot = FALSE)$acf[2:4]## [1] 0.1319241 -0.1870874 -0.2034868p=2 #orden del modelo
# sample autocorreleation function r
r=NULL
r[1:p]=acf(ar.process.LH, plot=F)$acf[2:(p+1)]
r## [1] 0.1319241 -0.1870874# matrix R
R=matrix(1,p,p) # matrix of dimension 4 by 4, with entries all 1's.
# define non-diagonal entires of R
for(i in 1:p){
for(j in 1:p){
if(i!=j)
R[i,j]=r[abs(i-j)]
}
}
cat('Matriz R de estimaciones de coeficientes de correlación.\n')## Matriz R de estimaciones de coeficientes de correlación.R## [,1] [,2]
## [1,] 1.0000000 0.1319241
## [2,] 0.1319241 1.0000000# b-column vector on the right
b=matrix(r,p,1)# b- column vector with no entries
cat('Vector b de estimaciones de coeficientes de correlación.\n')## Vector b de estimaciones de coeficientes de correlación.b## [,1]
## [1,] 0.1319241
## [2,] -0.1870874phi.hat=solve(R,b)[,1]
cat('Estimación de los coeficientes del modelo AR.\n')## Estimación de los coeficientes del modelo AR.phi.hat## [1] 0.1593793 -0.2081134# Variance estimation using Yule-Walker Estimator
c0=acf(ar.process.LH, type='covariance', plot=F)$acf[1]
cat('Estimación del coeficiente 0 de la autocovarianza.\n')## Estimación del coeficiente 0 de la autocovarianza.c0## [1] 0.5552905var.hat=c0*(1-sum(phi.hat*r))
cat('Valor de la varianza del ruido blando.\n')## Valor de la varianza del ruido blando.var.hat## [1] 0.5219945# Constant term in the model
phi0.hat=mean(df)*(1-sum(phi.hat))
cat('Estimación del coeficiente 0 del modelo AR.\n')## Estimación del coeficiente 0 del modelo AR.phi0.hat## [1] 607.2214cat('Modelo ajustado con las ecuaciones Yule-Walker para un proceso AR(',p,').\n')## Modelo ajustado con las ecuaciones Yule-Walker para un proceso AR( 2 ).cat("Constant:", phi0.hat," Coefficients:", phi.hat, " and Variance:", var.hat, '\n')## Constant: 607.2214 Coefficients: 0.1593793 -0.2081134 and Variance: 0.5219945Semana 5. Akaike Information Criterion (AIC), Mixed Models, Integrated Models
Akaike Information Criterion (AIC) and model quality
Vamos a simular con arima.sim() una serie para poder medir la calidad del modelo generado.
set.seed(43) #Roman conquest of Britain
data = arima.sim( list(order = c(2,0,0), ar =c( 0.7, -.2)), n = 2000) #AR(2) series
#hacemos las gracias del la serie, ACF y PACF
par(mfrow=c(1,3))
plot(data,main="AR(2) time series")
acf(data, main="ACF of AR Data of Second Order")
acf(data, type="partial", main="PACF of Time Series")
#el comando arima estima los coeficientes indicandole el orden del modelo
arima(data, order=c(2,0,0), include.mean=FALSE )##
## Call:
## arima(x = data, order = c(2, 0, 0), include.mean = FALSE)
##
## Coefficients:
## ar1 ar2
## 0.7111 -0.1912
## s.e. 0.0219 0.0220
##
## sigma^2 estimated as 0.9985: log likelihood = -2836.64, aic = 5679.27#también podemos ejecutarlo para varios ordenes
SSE=NULL
AIC=NULL
for (p in 1:5) {
m = arima(data, order=c(p,0,0), include.mean=FALSE )
SSE[p] = sum(resid(m)^2)
AIC[p] = m$aic
print(m$coef)
print(paste(m$aic, sum(resid(m)^2)))
}## ar1
## 0.5969948
## [1] "5751.73196762524 2072.83193501056"
## ar1 ar2
## 0.7111457 -0.1911552
## [1] "5679.27375222458 1997.00667996082"
## ar1 ar2 ar3
## 0.71359316 -0.20027407 0.01281966
## [1] "5680.94495534325 1996.67791506108"
## ar1 ar2 ar3 ar4
## 0.713676748 -0.201599647 0.017553046 -0.006629412
## [1] "5682.85704377107 1996.58997811231"
## ar1 ar2 ar3 ar4 ar5
## 0.71410825 -0.20268672 0.03019322 -0.05154692 0.06293048
## [1] "5676.91730818182 1988.65973372245"#Vamos a crear los graficos para el SSE y el AIC
par(mfrow=c(1,2))
order=c(1,2,3,4,5)
plot(SSE~order, main="SSE plotted on Order of AR(p) Model", ylim=c(1800, 2100))
plot(AIC~order, main="AIC plotted on Order of AR(p) Model", ylim=c(5500, 5800))
Siempre se prefiere un modelo con el AIC más bajo. Podemos ver en los graficos que el gran salto se da cuando pasamos de orden 1 a orden 2, consistente con un proceso AR(2).
ARMA (AutoRegressive Moving Average) models
En el mundo real, los procesos naturales tienen una parte autorregresiva (el valor actual tiene parte de los valores anteriores) y otra de medias moviles (el valor actual depende de las innovaciones/ruido blanco de los valores anteriores). PAra ello se usan los modelos mixtos como ARMA(p,q).
rm(list=ls(all=TRUE))
set.seed(500) # Beginning of Heptarchy: Kent, Essex, Sussex, Wessex, East Anglia, Mercia, and Northumbria.
data = arima.sim( list(order = c(1,0,1), ar =.7, ma=.2), n = 1000000) #ARMA(1,1) simulation 1M puntos
par(mfcol = c(3,1 ))
plot(data, main="ARMA(1,1) Time Series: phi1=.7, theta1=.2", xlim=c(0,4000)) #first 4000 terms
acf(data, main="Autocorrelation of ARMA(1,1), phi1=.7, theta1=.2")
acf(data, type="partial", main="Partial Autocorrelation of ARMA(1,1), phi1=.7, theta1=.2")
Para ver las propiedades de los procesos ARMA(p,q), vamos a usar la serie discoveries que lleva R por defecto.
plot(discoveries,
main = "Time Series of Number of Major Scientific Discoveries in a Year")
stripchart(discoveries, method = "stack", offset=.5, at=.15,pch=19,main="Number of Discoveries Dotplot",xlab="Number of Major Scientific Discoveries in a Year",ylab="Frequency")
#vamos a intentar averiguar el orden de los modelos ARMA con el acf y pacf
par(mfcol = c(2,1 ))
acf(discoveries, main="ACF of Number of Major Scientific Discoveries in a Year")
acf(discoveries, type="partial", main="PACF of Number of Major Scientific Discoveries
in a Year")
#A pesar que parece que el orden de MA es 1 o 2 o 3 (del ACF) y del AR es 1 o 2 también (del PACF), vamos a calcular el AIC de todas las combinaciones de p=0,1,2,3 y de q=0,1,2,3, y elegiremos el menor
AIC( arima( discoveries, order=c(0,0,1) ) ) #AIC = [1] 445.5895## [1] 445.5895AIC( arima( discoveries, order=c(0,0,2) ) ) #AIC = [1] 444.6742## [1] 444.6742AIC( arima( discoveries, order=c(0,0,3) ) ) #AIC = [1] 441.323## [1] 441.323AIC( arima( discoveries, order=c(1,0,0) ) ) #AIC = [1] 443.3792## [1] 443.3792AIC( arima( discoveries, order=c(1,0,1) ) ) #AIC = [1] 440.198## [1] 440.198AIC( arima( discoveries, order=c(1,0,2) ) ) #AIC = [1] 442.0428## [1] 442.0428AIC( arima( discoveries, order=c(1,0,3) ) ) #AIC = [1] 442.6747## [1] 442.6747AIC( arima( discoveries, order=c(2,0,0) ) ) #AIC = [1] 441.6155## [1] 441.6155AIC( arima( discoveries, order=c(2,0,1) ) ) #AIC = [1] 442.0722## [1] 442.0722AIC( arima( discoveries, order=c(2,0,2) ) ) #AIC = [1] 443.7021## [1] 443.7021AIC( arima( discoveries, order=c(2,0,3) ) ) #AIC = [1] 441.6594## [1] 441.6594AIC( arima( discoveries, order=c(3,0,0) ) ) #AIC = [1] 441.5658## [1] 441.5658AIC( arima( discoveries, order=c(3,0,1) ) ) #AIC = [1] 443.5655## [1] 443.5655AIC( arima( discoveries, order=c(3,0,2) ) ) #AIC = [1] 439.9263## [1] 439.9263AIC( arima( discoveries, order=c(3,0,3) ) ) #AIC = [1] 441.2941## [1] 441.2941#Los más bajos son el p=1 y q=1, y el p=3 y q=2. POr el principio de parsimonia, preferimos el ARMA(1,1), aunque el ARMA(3,2) tiene un AIC más bajo.
#vemos todos los detalles del modelo ARMA(1,1)
arima( discoveries, order=c(1,0,1) ) ##
## Call:
## arima(x = discoveries, order = c(1, 0, 1))
##
## Coefficients:
## ar1 ma1 intercept
## 0.8353 -0.6243 3.0208
## s.e. 0.1379 0.1948 0.4728
##
## sigma^2 estimated as 4.401: log likelihood = -216.1, aic = 440.2También podemos usar las rutinas automaticas del paquete forecast que vemos a continuación.
#cuidado con el flag approximation, se debe poner a TRUE si lenght(x) > 100 o frequency >12
#veamos como afecta a la estimación que nos da auto.arima
auto.arima(discoveries, d=0, approximation=FALSE)## Series: discoveries
## ARIMA(1,0,1) with non-zero mean
##
## Coefficients:
## ar1 ma1 mean
## 0.8353 -0.6243 3.0208
## s.e. 0.1379 0.1948 0.4728
##
## sigma^2 estimated as 4.538: log likelihood=-216.1
## AIC=440.2 AICc=440.62 BIC=450.62auto.arima(discoveries,d=0,approximation=TRUE)## Series: discoveries
## ARIMA(1,0,1) with non-zero mean
##
## Coefficients:
## ar1 ma1 mean
## 0.8353 -0.6243 3.0208
## s.e. 0.1379 0.1948 0.4728
##
## sigma^2 estimated as 4.538: log likelihood=-216.1
## AIC=440.2 AICc=440.62 BIC=450.62#tambien podemos elegir el criterio para elegir al modelo, entre Bayesian Information Criterion (BIC), Akaike Information Criterion (AIC), or “corrected AIC” (AICC). Por defecto está el Corrected AIC.
auto.arima(discoveries, d=0, ic="bic", approximation=FALSE)## Series: discoveries
## ARIMA(1,0,1) with non-zero mean
##
## Coefficients:
## ar1 ma1 mean
## 0.8353 -0.6243 3.0208
## s.e. 0.1379 0.1948 0.4728
##
## sigma^2 estimated as 4.538: log likelihood=-216.1
## AIC=440.2 AICc=440.62 BIC=450.62auto.arima(discoveries, d=0, ic="aic", approximation=FALSE)## Series: discoveries
## ARIMA(3,0,0) with non-zero mean
##
## Coefficients:
## ar1 ar2 ar3 mean
## 0.1967 0.1613 0.1451 3.0637
## s.e. 0.0995 0.0998 0.1007 0.4136
##
## sigma^2 estimated as 4.556: log likelihood=-215.78
## AIC=441.57 AICc=442.2 BIC=454.59auto.arima(discoveries, d=0, ic="aic", approximation=TRUE)## Series: discoveries
## ARIMA(1,0,1) with non-zero mean
##
## Coefficients:
## ar1 ma1 mean
## 0.8353 -0.6243 3.0208
## s.e. 0.1379 0.1948 0.4728
##
## sigma^2 estimated as 4.538: log likelihood=-216.1
## AIC=440.2 AICc=440.62 BIC=450.62ARIMA (AutoRegressive Integrated Moving Average) models
Hemos visto los modelos ARMA(p,q), donde hay una parte autorregresiva y una parte de media movil. Para estos modelos ARMA necesitamos que el proceso sea estacionario, para que la serie converja, pero es algo que no se da en la naturaleza. En la mayoria de los casos la no estacionariedad se debe a cierta tendencia en la serie. Debemos eliminar esa tendencia, y para ello usamos el operador diferencia, definido como delta=1-B, donde B es el backshift operator. La parte de tendencia es la que se introduce en la parte integrated de los modelos ARIMA.
Simulacion de un proceso ARIMA(2,1,1)
#parametros
phi=c(.7, .2)
beta=0.5
sigma=3
m=10000
#simulamos el proceso
set.seed(5)
Simulated.Arima <- arima.sim(n=m,list(order = c(2,1,1), ar = phi, ma=beta))
#hacemos el grafico de la serie
plot(Simulated.Arima, ylab=' ',main='Simulated time series from ARIMA(2,1,1) process', col='blue', lwd=2)
#vemos el ACF
acf(Simulated.Arima) #se ve que no es un proceso estacionario
#hacemos la diferencia en la serie, para eliminar la tendencia
Diff.Simulated.Arima <- diff(Simulated.Arima)
#vemos el grafico de la serie diferencia, el acf y el pacf
par(mfrow=c(3,1))
plot(Diff.Simulated.Arima)
acf(Diff.Simulated.Arima)
pacf(Diff.Simulated.Arima)
#buscamos el ajuste de un modelo a la serie original
sarima(Simulated.Arima,2,1,1,0,0,0) #fit ARIMA model## initial value 1.092704
## iter 2 value 0.655083
## iter 3 value 0.576329
## iter 4 value 0.250793
## iter 5 value 0.124855
## iter 6 value 0.033738
## iter 7 value 0.013225
## iter 8 value 0.012554
## iter 9 value 0.012517
## iter 10 value 0.012292
## iter 11 value 0.012267
## iter 12 value 0.012258
## iter 13 value 0.012170
## iter 14 value 0.012069
## iter 15 value 0.011860
## iter 16 value 0.011703
## iter 17 value 0.011609
## iter 18 value 0.011601
## iter 19 value 0.011601
## iter 20 value 0.011601
## iter 20 value 0.011601
## iter 20 value 0.011601
## final value 0.011601
## converged
## initial value 0.011653
## iter 2 value 0.011653
## iter 3 value 0.011653
## iter 3 value 0.011653
## iter 3 value 0.011653
## final value 0.011653
## converged
## $fit
##
## Call:
## stats::arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D,
## Q), period = S), xreg = constant, transform.pars = trans, fixed = fixed,
## optim.control = list(trace = trc, REPORT = 1, reltol = tol))
##
## Coefficients:
## ar1 ar2 ma1 constant
## 0.6876 0.204 0.5002 0.0280
## s.e. 0.0334 0.032 0.0301 0.1398
##
## sigma^2 estimated as 1.023: log likelihood = -14305.92, aic = 28621.83
##
## $degrees_of_freedom
## [1] 9996
##
## $ttable
## Estimate SE t.value p.value
## ar1 0.6876 0.0334 20.5786 0.0000
## ar2 0.2040 0.0320 6.3817 0.0000
## ma1 0.5002 0.0301 16.6139 0.0000
## constant 0.0280 0.1398 0.2001 0.8414
##
## $AIC
## [1] 2.862183
##
## $AICc
## [1] 2.862183
##
## $BIC
## [1] 2.865788auto.arima(Simulated.Arima) #evaluate best ARIMA model## Series: Simulated.Arima
## ARIMA(4,2,0)
##
## Coefficients:
## ar1 ar2 ar3 ar4
## 0.2279 -0.1633 0.0337 -0.0707
## s.e. 0.0100 0.0102 0.0102 0.0100
##
## sigma^2 estimated as 1.064: log likelihood=-14495.62
## AIC=29001.24 AICc=29001.25 BIC=29037.29#hacemos el ajuste a la serie diferencia, y vemos que pasa cuando no escogemos el modelo correcto
fit1 <- arima(Diff.Simulated.Arima, order=c(4,0,0))
fit1##
## Call:
## arima(x = Diff.Simulated.Arima, order = c(4, 0, 0))
##
## Coefficients:
## ar1 ar2 ar3 ar4 intercept
## 1.1862 -0.3761 0.1733 -0.0581 0.0280
## s.e. 0.0100 0.0154 0.0154 0.0100 0.1353
##
## sigma^2 estimated as 1.025: log likelihood = -14313.1, aic = 28638.2fit2 <- arima(Diff.Simulated.Arima, order=c(2,0,1))
fit2##
## Call:
## arima(x = Diff.Simulated.Arima, order = c(2, 0, 1))
##
## Coefficients:
## ar1 ar2 ma1 intercept
## 0.6876 0.204 0.5002 0.0280
## s.e. 0.0334 0.032 0.0301 0.1398
##
## sigma^2 estimated as 1.023: log likelihood = -14305.92, aic = 28621.83fit3 <- arima(Simulated.Arima, order=c(2,1,1)) #este es el orden correcto
fit3##
## Call:
## arima(x = Simulated.Arima, order = c(2, 1, 1))
##
## Coefficients:
## ar1 ar2 ma1
## 0.6876 0.2039 0.5001
## s.e. 0.0334 0.0320 0.0301
##
## sigma^2 estimated as 1.023: log likelihood = -14305.93, aic = 28619.85Ljung-Box Q-Statistics
Nos fijamos en el p-value de los residuos, que esperamos que sean ruido blando, y por tanto, los coeficientes de autocorrelación habn de ser cero. Esto ocurre con un valor p-value por debajo de un nivel de significación (tipicamente p-value < 0.05).
Ajuste de modelo a la serie temporal de nacimientos diarios
Vamos a hacer un ajuste de modelo a la serie temporal de los nacimientos diarios en california en 1959. TAmbien evaluaremos el mejor ajuste y usaremos herramientas conocidas como el ACF, PACF y AIC para discriminar entre modelos.
# The following time Series data is taken from Time Series Data Library (TSDL)
# TSDL was created by Rob Hyndman
# Professor of Statistics at Monash University, Australia.
# ==============================================================================
# ====== Daily total female birth in California, 1959 =======
# Data is exported as csv file to the wroking directory
# Link: https://datamarket.com/data/list/?q=cat:fwy%20provider:tsdl
par(mfrow=c(1,1))
rm(list=ls(all=TRUE))
# read data to R variable
birth.data <- read.csv("daily-total-female-births-in-cal.csv")
# pull out number of births column
number_of_births <- birth.data$Daily.total.female.births.in.California..1959
# use date format for dates
birth.data$Date <- as.Date(birth.data$Date, "%m/%d/%Y")
#dibujamos la serie
#tenemos tendencia y variación en la varianza
plot.ts(number_of_births,main='Daily total female births in california, 1959', ylab = 'Number of births')
# Comprobamos si hay correlación entre las muestras
Box.test(number_of_births, lag = log(length(number_of_births))) #al salir un p-value muy bajo, podemos rechazar la hipotesis nula de que no hay correlación, y por tanto, podemos extraer esta información para predecir valores futuros con un modelo como ARIMA##
## Box-Pierce test
##
## data: number_of_births
## X-squared = 36.391, df = 5.8999, p-value = 0.000002088# usamos la serie usando la diferencia, para intentar eliminar la tendencia en los datos, aunque tenemos picos, con la diferencia tenemos más o menos una serie estacionaria
plot.ts(diff(number_of_births), main='Differenced series', ylab = '')
# Con el test de autocorrelación nos sale un valor muy bajo de p-value, por lo que podemos rechazar la hipotesis nula y definitivamente hay correlación entre los datos que podemos extraer para hacer una predicción
Box.test(diff(number_of_births), lag = log(length(diff(number_of_births))))##
## Box-Pierce test
##
## data: diff(number_of_births)
## X-squared = 78.094, df = 5.8972, p-value = 0.000000000000007661# ACF y PACF de la serie diferencia, para poder estimar el orden de los procesos. Picos espureos en separaciones muy altas se pueden ignorar para el modelo, serán debidos a datos extremos o aleatorios
acf(diff(number_of_births), main='ACF of differenced data', 50) #este ACF sugiere un proceso MA(1)
pacf(diff(number_of_births), main='PACF of differenced data', 50) #este PACF sugiere un proceso AR(7)
# Ajustamos varios modelos ARIMA a la serie original. Elegimos d=1 en todos porque hemos hecho una diferencia en la serie para estimar los ordenes de los procesos AR y MA involucrados
model1 <- arima(number_of_births, order=c(0,1,1))
SSE1 <- sum(model1$residuals^2)
model1.test <- Box.test(model1$residuals, lag = log(length(model1$residuals)))
model2 <- arima(number_of_births, order=c(0,1,2))
SSE2 <- sum(model2$residuals^2)
model2.test <- Box.test(model2$residuals, lag = log(length(model2$residuals)))
model3 <- arima(number_of_births, order=c(7,1,1))
SSE3 <- sum(model3$residuals^2)
model3.test <- Box.test(model3$residuals, lag = log(length(model3$residuals)))
model4 <- arima(number_of_births, order=c(7,1,2))
SSE4 <- sum(model4$residuals^2)
model4.test <- Box.test(model4$residuals, lag = log(length(model4$residuals)))
df <- data.frame(row.names=c('AIC', 'SSE', 'p-value'), c(model1$aic, SSE1, model1.test$p.value),
c(model2$aic, SSE2, model2.test$p.value), c(model3$aic, SSE3, model3.test$p.value),
c(model4$aic, SSE4, model4.test$p.value))
colnames(df) <- c('Arima(0,1,1)','Arima(0,1,2)', 'Arima(7,1,1)', 'Arima(7,1,2)')
format(df, scientific=FALSE)## Arima(0,1,1) Arima(0,1,2) Arima(7,1,1) Arima(7,1,2)
## AIC 2462.2207021 2459.5705306 2464.8827225 2466.6664136
## SSE 18148.4561632 17914.6513437 17584.3902548 17574.0578107
## p-value 0.5333604 0.9859227 0.9999899 0.9999929#En este punto, como todos los modelos tienen los residuos sin autocorrelacion (son ruido blanco), debemos elegir el modelo con menor AIC o SSE. El menor AIC lo tiene el modelo ARIMA(0,1,2) y el menor SSE lo tiene el modelo ARIMA(7,1,2). Este ultimo tendría 10 parametros y seria complejo, siempre elegiremos el modelo más simple, por lo que nos inclinaremos por el AIC menos y el modelo ARIMA(0,1,2) para modelar estos datos.
# Usamos SARIMA para estimar el mocelo y los coeficientes y otros estaditicos con el modelo elegido
sarima(number_of_births,0,1,2,0,0,0)## initial value 2.216721
## iter 2 value 2.047518
## iter 3 value 1.974780
## iter 4 value 1.966955
## iter 5 value 1.958906
## iter 6 value 1.952299
## iter 7 value 1.951439
## iter 8 value 1.950801
## iter 9 value 1.950797
## iter 10 value 1.950650
## iter 11 value 1.950646
## iter 12 value 1.950638
## iter 13 value 1.950635
## iter 13 value 1.950635
## iter 13 value 1.950635
## final value 1.950635
## converged
## initial value 1.950708
## iter 2 value 1.950564
## iter 3 value 1.950290
## iter 4 value 1.950196
## iter 5 value 1.950185
## iter 6 value 1.950185
## iter 7 value 1.950185
## iter 7 value 1.950185
## iter 7 value 1.950185
## final value 1.950185
## converged
## $fit
##
## Call:
## stats::arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D,
## Q), period = S), xreg = constant, transform.pars = trans, fixed = fixed,
## optim.control = list(trace = trc, REPORT = 1, reltol = tol))
##
## Coefficients:
## ma1 ma2 constant
## -0.8511 -0.1113 0.015
## s.e. 0.0496 0.0502 0.015
##
## sigma^2 estimated as 49.08: log likelihood = -1226.36, aic = 2460.72
##
## $degrees_of_freedom
## [1] 361
##
## $ttable
## Estimate SE t.value p.value
## ma1 -0.8511 0.0496 -17.1448 0.0000
## ma2 -0.1113 0.0502 -2.2164 0.0273
## constant 0.0150 0.0150 1.0007 0.3176
##
## $AIC
## [1] 6.760225
##
## $AICc
## [1] 6.760408
##
## $BIC
## [1] 6.803051Por lo tanto, con la salida de sarima podemos estimar nuestro proceso como:
(1-B)Xt = 0.015 + Zt - 0.8511 * Z(t-1) -0.1113 * Z(t-2)
y operando el operador backshift:
Xt=X(t-1) + Zt - 0.8511 * Z(t-1) -0.1113 * Z(t-2)
donde Zt es un ruido blanco de media 0 y varianza de 49.08. Con esta expresion matematica ya podriamos hacer una predicción de esta serie temporal.
Ajuste de la serie BJsales contenida en R
# Plot time series 'BJsales'
plot(BJsales,main="Time series BJsales")
plot(diff(BJsales)) #continua habiendo tendencia en los datos
plot(diff(diff(BJsales))) #tendremos un modelo ARIMA con d=2
acf(diff(diff(BJsales)),type="partial") #PACF, sugiere AR de orden 0,1,2,3
acf(diff(diff(BJsales))) #ACF, sugiere MA de orden 0 o 1
#probamos el ajuste de varios modelos ARIMA, el de menor AIC es el ARIMA(0,2,1) y el de menor SSE es el ARIMA(3,2,1)
d=2
for(p in 1:4){
for(q in 1:2){
if(p+d+q<=8){
model <- arima(x=BJsales, order = c((p-1),d,(q-1)))
pval <- Box.test(model$residuals, lag=log(length(model$residuals)))
sse <- sum(model$residuals^2)
cat(p-1,d,q-1, 'AIC=', model$aic, ' SSE=',sse,' p-VALUE=', pval$p.value,'\n')
}
}
}## 0 2 0 AIC= 577.6777 SSE= 423.7908 p-VALUE= 0.0000007610494
## 0 2 1 AIC= 517.1371 SSE= 276.2293 p-VALUE= 0.9632467
## 1 2 0 AIC= 541.9646 SSE= 327.92 p-VALUE= 0.003606983
## 1 2 1 AIC= 518.9734 SSE= 275.8554 p-VALUE= 0.941776
## 2 2 0 AIC= 532.2986 SSE= 302.7467 p-VALUE= 0.05824465
## 2 2 1 AIC= 520.2684 SSE= 274.0474 p-VALUE= 0.7955439
## 3 2 0 AIC= 524.7648 SSE= 283.4941 p-VALUE= 0.7035291
## 3 2 1 AIC= 519.4182 SSE= 264.0684 p-VALUE= 0.6948066#por el principio de parsimonia, elegimos el de menor AIC, el ARIMA(0,2,1)
#evaluamos los residuos de este modelo, para comprobar que son ruido blanco
model <- arima(BJsales, order=c(0,2,1))
par(mfrow=c(2,2))
plot(model$residuals)
acf(model$residuals)
pacf(model$residuals)
qqnorm(model$residuals)
#obtenemos los coeficientes del modelo con sarima
sarima(BJsales,0,2,1,0,0,0)## initial value 0.525918
## iter 2 value 0.353629
## iter 3 value 0.330007
## iter 4 value 0.329249
## iter 5 value 0.315928
## iter 6 value 0.313389
## iter 7 value 0.312977
## iter 8 value 0.312970
## iter 9 value 0.312965
## iter 9 value 0.312965
## iter 9 value 0.312965
## final value 0.312965
## converged
## initial value 0.314633
## iter 1 value 0.314633
## final value 0.314633
## converged
## $fit
##
## Call:
## stats::arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D,
## Q), period = S), include.mean = !no.constant, transform.pars = trans, fixed = fixed,
## optim.control = list(trace = trc, REPORT = 1, reltol = tol))
##
## Coefficients:
## ma1
## -0.7480
## s.e. 0.0662
##
## sigma^2 estimated as 1.866: log likelihood = -256.57, aic = 517.14
##
## $degrees_of_freedom
## [1] 147
##
## $ttable
## Estimate SE t.value p.value
## ma1 -0.748 0.0662 -11.3045 0
##
## $AIC
## [1] 3.49417
##
## $AICc
## [1] 3.494355
##
## $BIC
## [1] 3.534672Por lo tanto, el modelo será el siguiente:
(1-B)^2 Xt = (1-0.748) * Zt
Xt = 2 X(t-1) + X(t-2) + Zt -0.748 Z(t-2)
donde Zt es un ruido blando de media 0 y varianza 1.866.
Semana 6 - Estacionalidad, SARIMA y Predicción (Forecasting)
SARIMA
En series reales, a los componentes autorregresivos y de medias moviles se le añade un componente de estacionalidad, componentes que aparecen cada cierto tiempo, que son periódicos cada cierto numero de observaciones (por ejemplo, en series mensuales s=12, en series trimestrales, s=4). Para modelar estas series tenemos los modelos SARIMA (Seasonal ARIMA).
Los modelos SARIMA tienen 7 variables - SARIMA (p,d,q,P,D,Q)s
- p - orden de la componente autorregresiva, AR(p)
- d - orden de la diferencia
- q - orden de la componente de medias moviles, MA(q)
- P - orden de la componente autorregresiva estacional, SAR(P)
- D - orden de la diferencia estacional
- Q - orden de la componente de medias moviles estacional, SMA(Q)
- s - parámetro de estacionalidad
VAmos a simular un modelo SARIMA(0,0,1,0,0,1)12, y vemos como se comporta su ACF.
x=NULL
z=NULL
n=10000
z=rnorm(n)
x[1:13]=1
for(i in 14:n){
x[i] <- z[i]+0.7*z[i-1]+0.6*z[i-12]+0.42*z[i-13]
}
plot.ts(x[1:2000],main="SARIMA(0,0,1,0,0,1)12 model - 2000 points")
par(mfrow=c(2,1))
plot.ts(x[12:120], main='The first 10 months of simulation SARIMA(0,0,1,0,0)_12', ylab='')
acf(x,main='SARIMA(0,0,1,0,0,1)_12 Simulation')
En el ACF vemos el pico en la diferencia 1, lo que es normal en un proceso con componente MA(1), y además vemos los picos alrededor del parametro de estacionalidad, en este caso s=12 y en 13, tal y como explica el modelo matemático. Pero también tenemos un pico en 11 en el ACF, que no aparece directamente en el modelo por ninguna parte. Esto sucede porque la autocorrelacion en 11 dependerá también de Z12, y por tanto, no será cero o poco significativo.
Aplicaciones de los modelos SARIMA
Pasos para el modelado de una serie temporal:
- Dibujar gráficamente la serie.
- Transformación de la serie (la más común son los logaritmos) para hacerla mas estacional, sobre todo en varianza.
- Diferencia (estacional o no estacional) para eliminar la tendencia - nos da los valores d y D.
- Ljung-Box test, para ver si hay alguna correlación entre muestras y por tanto, podemos modelar la serie (p-value < 0.05).
- ACF - picos significativos - orden q del proceso MA(q).
- ACF - picos significativos en alrededor de periodos estacionales - orden Q del proceso SMA(Q).
- PACF - picos significativos - orden p del proceso AR(p).
- PACF - picos significativos en alrededor de periodos estacionales - orden P del proceso SAR(P).
- Ajuste con ARIMA para la elección del modelo con menor AIC o SSE, y usando el principio de parsimonia (p+d+q+P+D+Q<=6). Este principio evita modelos muy complicados que pueden sobreestimar el modelo.
- Modelar con SARIMA con el modelo elegido para calcular la estimación de los coeficientes.
- Dibujar graficamente los residuos, su ACF y su PACF y ver que no hay valores significativos.
- Comprobación de la aleatoriedad de los residuos con Ljung-Box test (p-value > 0.05).
- Usar forecast para estimar valores futuros según el modelo.
Ajuste de la serie de J&J con SARIMA
#inicializamos parametros
d=1
DD=1
per=4
#modelado con ARIMA para elegir el mejor modelo con menor AIC
for(p in 1:2){
for(q in 1:2){
for(i in 1:2){
for(j in 1:2){
if(p+d+q+i+DD+j<=10){
model <- arima(x=log(jj), order = c((p-1),d,(q-1)), seasonal = list(order=c((i-1),DD,(j-1)), period=per))
pval <- Box.test(model$residuals, lag=log(length(model$residuals)))
sse <- sum(model$residuals^2)
cat(p-1,d,q-1,i-1,DD,j-1,per, 'AIC=', model$aic, ' SSE=',sse,' p-VALUE=', pval$p.value,'\n')
}
}
}
}
}## 0 1 0 0 1 0 4 AIC= -124.0685 SSE= 0.9377871 p-VALUE= 0.0002610795
## 0 1 0 0 1 1 4 AIC= -126.3493 SSE= 0.8856994 p-VALUE= 0.0001606503
## 0 1 0 1 1 0 4 AIC= -125.9198 SSE= 0.8908544 p-VALUE= 0.0001978064
## 0 1 0 1 1 1 4 AIC= -124.3648 SSE= 0.8854554 p-VALUE= 0.000157403
## 0 1 1 0 1 0 4 AIC= -145.5139 SSE= 0.6891988 p-VALUE= 0.03543716
## 0 1 1 0 1 1 4 AIC= -150.7528 SSE= 0.6265214 p-VALUE= 0.6089542
## 0 1 1 1 1 0 4 AIC= -150.9134 SSE= 0.6251634 p-VALUE= 0.7079181
## 0 1 1 1 1 1 4 AIC= -149.1317 SSE= 0.6232876 p-VALUE= 0.6780876
## 1 1 0 0 1 0 4 AIC= -139.8248 SSE= 0.7467494 p-VALUE= 0.03503415
## 1 1 0 0 1 1 4 AIC= -146.0191 SSE= 0.6692691 p-VALUE= 0.5400211
## 1 1 0 1 1 0 4 AIC= -146.0319 SSE= 0.6689661 p-VALUE= 0.5612964
## 1 1 0 1 1 1 4 AIC= -144.3766 SSE= 0.6658382 p-VALUE= 0.5459445
## 1 1 1 0 1 0 4 AIC= -145.8284 SSE= 0.667109 p-VALUE= 0.2200484
## 1 1 1 0 1 1 4 AIC= -148.7706 SSE= 0.6263677 p-VALUE= 0.5948219
## 1 1 1 1 1 0 4 AIC= -148.9175 SSE= 0.6251104 p-VALUE= 0.719547
## 1 1 1 1 1 1 4 AIC= -144.4483 SSE= 0.6097742 p-VALUE= 0.3002703#modelado con SARIMA
sarima(log(jj), 0,1,1,1,1,0,4)## initial value -2.237259
## iter 2 value -2.429075
## iter 3 value -2.446738
## iter 4 value -2.455821
## iter 5 value -2.459761
## iter 6 value -2.462511
## iter 7 value -2.462602
## iter 8 value -2.462749
## iter 9 value -2.462749
## iter 9 value -2.462749
## iter 9 value -2.462749
## final value -2.462749
## converged
## initial value -2.411490
## iter 2 value -2.412022
## iter 3 value -2.412060
## iter 4 value -2.412062
## iter 4 value -2.412062
## iter 4 value -2.412062
## final value -2.412062
## converged
## $fit
##
## Call:
## stats::arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D,
## Q), period = S), include.mean = !no.constant, transform.pars = trans, fixed = fixed,
## optim.control = list(trace = trc, REPORT = 1, reltol = tol))
##
## Coefficients:
## ma1 sar1
## -0.6796 -0.3220
## s.e. 0.0969 0.1124
##
## sigma^2 estimated as 0.007913: log likelihood = 78.46, aic = -150.91
##
## $degrees_of_freedom
## [1] 77
##
## $ttable
## Estimate SE t.value p.value
## ma1 -0.6796 0.0969 -7.0104 0.0000
## sar1 -0.3220 0.1124 -2.8641 0.0054
##
## $AIC
## [1] -1.840408
##
## $AICc
## [1] -1.838555
##
## $BIC
## [1] -1.753721#dibujamos la predicción del modelo elegido
model <- arima(x=log(jj),order=c(0,1,1),seasonal=list(order=c(1,1,0),period=4))
par(mfrow=c(1,1))
plot(forecast(model))
forecast(model) #estimación de los siguientes 8 valores (2 años)## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 1981 Q1 2.910254 2.796250 3.024258 2.735900 3.084608
## 1981 Q2 2.817218 2.697507 2.936929 2.634135 3.000300
## 1981 Q3 2.920738 2.795579 3.045896 2.729325 3.112151
## 1981 Q4 2.574797 2.444419 2.705175 2.375401 2.774194
## 1982 Q1 3.041247 2.868176 3.214317 2.776558 3.305935
## 1982 Q2 2.946224 2.762623 3.129824 2.665431 3.227016
## 1982 Q3 3.044757 2.851198 3.238316 2.748734 3.340780
## 1982 Q4 2.706534 2.503505 2.909564 2.396028 3.017041Ajuste de la serie temporal de producción de leche
milk <- read.csv('monthly-milk-production-pounds.csv')
names(milk)[2] <- 'Pounds'
Milk <- milk$Pounds
sarima(Milk, 0,1,0,0,1,1,12)## initial value 2.187412
## iter 2 value 2.041962
## iter 3 value 2.026713
## iter 4 value 2.023363
## iter 5 value 2.020977
## iter 6 value 2.020123
## iter 7 value 2.020102
## iter 7 value 2.020102
## final value 2.020102
## converged
## initial value 2.027140
## iter 2 value 2.027140
## iter 3 value 2.027140
## iter 3 value 2.027140
## iter 3 value 2.027140
## final value 2.027140
## converged
## $fit
##
## Call:
## stats::arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D,
## Q), period = S), include.mean = !no.constant, transform.pars = trans, fixed = fixed,
## optim.control = list(trace = trc, REPORT = 1, reltol = tol))
##
## Coefficients:
## sma1
## -0.6110
## s.e. 0.0622
##
## sigma^2 estimated as 55.6: log likelihood = -534.14, aic = 1072.28
##
## $degrees_of_freedom
## [1] 154
##
## $ttable
## Estimate SE t.value p.value
## sma1 -0.611 0.0622 -9.8265 0
##
## $AIC
## [1] 6.459545
##
## $AICc
## [1] 6.459692
##
## $BIC
## [1] 6.496212d=NULL
DD=NULL
d=1
DD=1
per=12
for(p in 1:1){
for(q in 1:1){
for(i in 1:3){
for(j in 1:4){
if(p+d+q+i+DD+j<=10){
model<-arima(x=Milk, order = c((p-1),d,(q-1)), seasonal = list(order=c((i-1),DD,(j-1)), period=per))
pval<-Box.test(model$residuals, lag=log(length(model$residuals)))
sse<-sum(model$residuals^2)
cat(p-1,d,q-1,i-1,DD,j-1,per, 'AIC=', model$aic, ' SSE=',sse,' p-VALUE=', pval$p.value,'\n')
}
}
}
}
}## 0 1 0 0 1 0 12 AIC= 1119.969 SSE= 12315.74 p-VALUE= 0.02571055
## 0 1 0 0 1 1 12 AIC= 1072.284 SSE= 8622.065 p-VALUE= 0.02663598
## 0 1 0 0 1 2 12 AIC= 1074.091 SSE= 8604.425 p-VALUE= 0.02484041
## 0 1 0 0 1 3 12 AIC= 1074.142 SSE= 8433.924 p-VALUE= 0.04199235
## 0 1 0 1 1 0 12 AIC= 1089.116 SSE= 9804.661 p-VALUE= 0.01144098
## 0 1 0 1 1 1 12 AIC= 1074.136 SSE= 8608.948 p-VALUE= 0.02510353
## 0 1 0 1 1 2 12 AIC= 1075.003 SSE= 8418.505 p-VALUE= 0.02820141
## 0 1 0 1 1 3 12 AIC= 1076.958 SSE= 8405.012 p-VALUE= 0.02994421
## 0 1 0 2 1 0 12 AIC= 1080.89 SSE= 9081.87 p-VALUE= 0.02797638
## 0 1 0 2 1 1 12 AIC= 1075.405 SSE= 8553.762 p-VALUE= 0.03308289
## 0 1 0 2 1 2 12 AIC= 1076.995 SSE= 8442.979 p-VALUE= 0.02902194#ajustamos el modelo
model<- arima(x=Milk, order = c(0,1,0), seasonal = list(order=c(0,1,1), period=12))
#hacemos la predicción de dos años segun el ajuste del modelo
plot(forecast(model))
forecast(model)## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 169 866.6448 857.0892 876.2003 852.0308 881.2588
## 170 819.2455 805.7319 832.7591 798.5782 839.9128
## 171 926.0938 909.5431 942.6446 900.7817 951.4060
## 172 939.1004 919.9892 958.2115 909.8724 968.3283
## 173 1002.2647 980.8978 1023.6316 969.5868 1034.9426
## 174 974.8489 951.4426 998.2552 939.0521 1010.6457
## 175 933.5422 908.2605 958.8239 894.8772 972.2072
## 176 893.9598 866.9325 920.9870 852.6252 935.2944
## 177 848.0580 819.3913 876.7247 804.2160 891.9000
## 178 853.2545 823.0372 883.4719 807.0410 899.4680
## 179 819.2343 787.5421 850.9266 770.7652 867.7034
## 180 861.5723 828.4708 894.6738 810.9479 912.1966
## 181 885.2171 849.5536 920.8805 830.6746 939.7595
## 182 837.8178 799.7645 875.8710 779.6203 896.0152
## 183 944.6661 904.3645 984.9678 883.0301 1006.3022
## 184 957.6726 915.2416 1000.1037 892.7800 1022.5653
## 185 1020.8370 976.3784 1065.2956 952.8435 1088.8305
## 186 993.4212 947.0236 1039.8188 922.4622 1064.3802
## 187 952.1145 903.8557 1000.3733 878.3091 1025.9199
## 188 912.5321 862.4813 962.5828 835.9860 989.0781
## 189 866.6303 814.8495 918.4111 787.4384 945.8222
## 190 871.8268 818.3720 925.2817 790.0747 953.5790
## 191 837.8066 782.7285 892.8847 753.5720 922.0412
## 192 880.1446 823.4898 936.7994 793.4985 966.7906Ajuste de la serie temporal de ventas de una tienda de souvenirs
SUV<-read.csv('fancy.csv')
suv<-ts(SUV$Sales)
#dibujamos graficamente las diferentes series temporales
par(mfrow=c(2,2))
plot(suv, main='Monthly sales for a souvenir shop', ylab='', col='blue', lwd=3)
plot(log(suv), main='Log-transorm of sales', ylab='', col='red', lwd=3)
plot(diff(log(suv)), main='Differenced Log-transorm of sales', ylab='', col='brown', lwd=3)
plot(diff(diff(log(suv)),12), main='Log-transorm without trend and seasonaliy', ylab='', col='green', lwd=3)
#transformamos los datos
data <- diff(diff((log(suv)),12))
acf2(data, 50)
## ACF PACF
## [1,] -0.46 -0.46
## [2,] 0.19 -0.02
## [3,] -0.17 -0.11
## [4,] -0.06 -0.23
## [5,] 0.01 -0.13
## [6,] 0.00 -0.07
## [7,] -0.07 -0.20
## [8,] 0.07 -0.12
## [9,] 0.09 0.11
## [10,] 0.02 0.11
## [11,] 0.00 0.04
## [12,] -0.20 -0.20
## [13,] 0.08 -0.07
## [14,] -0.05 -0.01
## [15,] 0.07 -0.01
## [16,] -0.03 -0.05
## [17,] 0.05 0.02
## [18,] -0.02 -0.03
## [19,] -0.05 -0.19
## [20,] 0.15 0.13
## [21,] -0.24 -0.06
## [22,] 0.18 0.01
## [23,] -0.01 0.15
## [24,] -0.07 -0.09
## [25,] 0.13 0.03
## [26,] -0.21 -0.16
## [27,] 0.15 0.04
## [28,] -0.11 -0.05
## [29,] 0.12 0.06
## [30,] -0.01 0.10
## [31,] -0.05 -0.12
## [32,] -0.04 -0.12
## [33,] 0.14 0.08
## [34,] -0.25 -0.15
## [35,] 0.15 -0.08
## [36,] -0.10 -0.05
## [37,] 0.09 0.00
## [38,] 0.05 -0.07
## [39,] 0.01 0.00
## [40,] -0.05 -0.07
## [41,] -0.04 0.01
## [42,] -0.07 -0.21
## [43,] 0.10 0.03
## [44,] -0.01 0.13
## [45,] 0.03 -0.05
## [46,] 0.02 -0.01
## [47,] 0.03 0.02
## [48,] -0.03 0.04
## [49,] 0.00 -0.02
## [50,] -0.05 -0.02#buscamos el modelo que más se ajuste
d=1
DD=1
per=12
for(p in 1:2){
for(q in 1:2){
for(i in 1:2){
for(j in 1:4){
if(p+d+q+i+DD+j<=10){
model<-arima(x=log(suv), order = c((p-1),d,(q-1)), seasonal = list(order=c((i-1),DD,(j-1)), period=per))
pval<-Box.test(model$residuals, lag=log(length(model$residuals)))
sse<-sum(model$residuals^2)
cat(p-1,d,q-1,i-1,DD,j-1,per, 'AIC=', model$aic, ' SSE=',sse,' p-VALUE=', pval$p.value,'\n')
}
}
}
}
}## 0 1 0 0 1 0 12 AIC= -11.60664 SSE= 3.432906 p-VALUE= 0.0001365568
## 0 1 0 0 1 1 12 AIC= -16.09179 SSE= 2.977559 p-VALUE= 0.00003149961
## 0 1 0 0 1 2 12 AIC= -17.58234 SSE= 2.301961 p-VALUE= 0.0002456592
## 0 1 0 0 1 3 12 AIC= -16.41016 SSE= 2.35266 p-VALUE= 0.0003392286
## 0 1 0 1 1 0 12 AIC= -13.43083 SSE= 3.214065 p-VALUE= 0.00004083829
## 0 1 0 1 1 1 12 AIC= -17.76362 SSE= 2.399746 p-VALUE= 0.0001916567
## 0 1 0 1 1 2 12 AIC= -15.99095 SSE= 2.349899 p-VALUE= 0.0002477785
## 0 1 0 1 1 3 12 AIC= -14.74777 SSE= 2.302026 p-VALUE= 0.0004504603
## 0 1 1 0 1 0 12 AIC= -27.78538 SSE= 2.643277 p-VALUE= 0.1742481
## 0 1 1 0 1 1 12 AIC= -34.54538 SSE= 2.233424 p-VALUE= 0.2730773
## 0 1 1 0 1 2 12 AIC= -33.6145 SSE= 2.10948 p-VALUE= 0.2830601
## 0 1 1 0 1 3 12 AIC= -32.19273 SSE= 1.87789 p-VALUE= 0.2700422
## 0 1 1 1 1 0 12 AIC= -32.33191 SSE= 2.360508 p-VALUE= 0.2584529
## 0 1 1 1 1 1 12 AIC= -34.0881 SSE= 1.842013 p-VALUE= 0.2843227
## 0 1 1 1 1 2 12 AIC= -32.1017 SSE= 1.856344 p-VALUE= 0.2851601
## 1 1 0 0 1 0 12 AIC= -27.07825 SSE= 2.6747 p-VALUE= 0.2297873
## 1 1 0 0 1 1 12 AIC= -34.98918 SSE= 2.209442 p-VALUE= 0.4633773
## 1 1 0 0 1 2 12 AIC= -33.38623 SSE= 2.159413 p-VALUE= 0.4515446
## 1 1 0 0 1 3 12 AIC= -31.54519 SSE= 2.121635 p-VALUE= 0.4390831
## 1 1 0 1 1 0 12 AIC= -32.64858 SSE= 2.340077 p-VALUE= 0.4022225
## 1 1 0 1 1 1 12 AIC= -33.48894 SSE= 2.125764 p-VALUE= 0.4442667
## 1 1 0 1 1 2 12 AIC= -31.52137 SSE= 2.093127 p-VALUE= 0.4463098
## 1 1 1 0 1 0 12 AIC= -26.17089 SSE= 2.624282 p-VALUE= 0.2507445
## 1 1 1 0 1 1 12 AIC= -33.30647 SSE= 2.201798 p-VALUE= 0.4110141
## 1 1 1 0 1 2 12 AIC= -31.68924 SSE= 2.151774 p-VALUE= 0.3820817
## 1 1 1 1 1 0 12 AIC= -31.10127 SSE= 2.323818 p-VALUE= 0.3492748
## 1 1 1 1 1 1 12 AIC= -32.69913 SSE= 1.823509 p-VALUE= 0.3092363model<- arima(x=log(suv), order = c(1,1,0), seasonal = list(order=c(0,1,1), period=12))
par(mfrow=c(1,1))
plot(forecast(model))
forecast(model)## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 85 9.600019 9.373967 9.826071 9.254303 9.945736
## 86 9.786505 9.533944 10.039066 9.400246 10.172763
## 87 10.329604 10.025423 10.633786 9.864399 10.794810
## 88 10.081972 9.746705 10.417239 9.569226 10.594719
## 89 10.008096 9.638604 10.377587 9.443008 10.573184
## 90 10.181170 9.783094 10.579245 9.572366 10.789973
## 91 10.439372 10.013362 10.865383 9.787846 11.090899
## 92 10.534857 10.083238 10.986477 9.844165 11.225550
## 93 10.613026 10.136886 11.089165 9.884834 11.341218
## 94 10.664526 10.165207 11.163845 9.900884 11.428169
## 95 11.096784 10.575248 11.618320 10.299164 11.894405
## 96 11.877167 11.334355 12.419978 11.047008 12.707325
## 97 9.932756 9.330374 10.535138 9.011492 10.854020
## 98 10.112193 9.475681 10.748706 9.138732 11.085655
## 99 10.658829 9.980845 11.336813 9.621941 11.695716
## 100 10.409423 9.696789 11.122057 9.319543 11.499303
## 101 10.336436 9.588667 11.084205 9.192822 11.480051
## 102 10.509064 9.728750 11.289377 9.315677 11.702450
## 103 10.767491 9.955450 11.579531 9.525582 12.009399
## 104 10.862863 10.020525 11.705201 9.574618 12.151108
## 105 10.941088 10.069391 11.812785 9.607942 12.274234
## 106 10.992560 10.092517 11.892604 9.616063 12.369058
## 107 11.424832 10.497281 12.352384 10.006265 12.843400
## 108 12.205208 11.250956 13.159460 10.745805 13.664611a <- sarima.for(log(suv),12,1,1,0,0,1,1,12)
plot.ts(c(suv,exp(a$pred)), main='Monthly sales + Forecast', ylab='', col='blue', lwd=3)
Forecasting con técnicas de suavizado
Vamos a explorar otras formas de estimar los valores futuros basandose en los valores pasados de una serie temporal, aparte de la vista hasta ahora. ¿Que podemos hacer, por ejemplo, si los modelos ARIMA no se ajustan bien a la serie temporal que tenemos?
rm(list=ls(all=TRUE))
rain.data <- scan("http://robjhyndman.com/tsdldata/hurst/precip1.dat",skip=1)
write.csv(rain.data,file = "rain_data.csv")
rain.ts <- ts(rain.data, start=c(1813))
#algunos graficos de la serie, vemos que no hay una distribución normal
par( mfrow=c(1,2) )
hist(rain.data, main="Annual London Rainfall 1813-1912",
xlab="rainfall in inches")
qqnorm(rain.data,main="Normal Plot of London Rainfall")
qqline(rain.data)
#dibujamos la serie y el ACF
par( mfrow=c(2,1) )
plot.ts(rain.ts, main="Annual London Rainfall 1813-1912",xlab="year", ylab="rainfall in inches")
acf(rain.ts, main="ACF: London Rainfall") #no parece que haya unas correlaciones fuertes como para hacer un modelo, a pesar que no es una distribución normal
#probamos con autoarima
auto.arima(rain.ts) #nos da un modelo 0,0,0 porque en realidad no encuentra ningún modelo## Series: rain.ts
## ARIMA(0,0,0) with non-zero mean
##
## Coefficients:
## mean
## 24.8239
## s.e. 0.4193
##
## sigma^2 estimated as 17.76: log likelihood=-285.25
## AIC=574.49 AICc=574.61 BIC=579.7Suavizado exponencial simple - SES
El suavizado exponencial simple (SES - Simple Exponential Smoothing) es una forma de estimar el siguiente valor dando más peso a los valores recientes, y asi se mejora un metodo de simplemente estimar la media de unos cuantos valores anteriores. Los pesos serán una serie geometrica convergente, donde los valores decaen. Elegir la razon de la serie geomética se ha de hacer en base a los errores que cometemos, y elegir aquella razón que los minimize.
#simulamos el SES con un bucle for
alpha <- .2 #increase alpha for more rapid decay
forecast.values <- NULL #establish array to store forecast values
n <- length(rain.data)
#naive first forecast
forecast.values [1] <- rain.data[1]
#loop to create all forecast values
for( i in 1:n ) {
forecast.values [i+1] <- alpha*rain.data[i] + (1-alpha)* forecast.values [i]
}
paste("forecast for time",n+1," = ", forecast.values [n+1])## [1] "forecast for time 101 = 25.3094062064236"#calculamos el SSE para varios valores de la razon alpha de la serie geométrica, para elegir el optimom que minimice los errores
SSE <- NULL
n <- length(rain.data)
alpha.values <- seq( .001, .999, by=0.001)
number.alphas <- length(alpha.values)
for( k in 1:number.alphas ) {
forecast.values <-NULL
alpha <- alpha.values[k]
forecast.values[1] <- rain.data[1]
for( i in 1:n ) {
forecast.values[i+1] <- alpha*rain.data[i] + (1-alpha)*forecast.values[i]
}
SSE[k] <- sum( (rain.data - forecast.values[1:n])^2 )
}
plot(SSE~alpha.values, main="Optimal alpha value Minimizes SSE")
index.of.smallest.SSE <- which.min(SSE) #returns position 24
alpha.values[which.min(SSE)] #este es el alpha que minimiza el error## [1] 0.024#hacemos la predicción con el valor de alpha optimo
alpha <- alpha.values[which.min(SSE)] #increase alpha for more rapid decay
forecast.values <- NULL #establish array to store forecast values
n <- length(rain.data)
#naive first forecast
forecast.values [1] <- rain.data[1]
#loop to create all forecast values
for( i in 1:n ) {
forecast.values [i+1] <- alpha*rain.data[i] + (1-alpha)* forecast.values [i]
}
paste("forecast for time",n+1," = ", forecast.values [n+1])## [1] "forecast for time 101 = 24.6771392918524"#comparamos nuestro calculo con la funcion Holt-Winters que veremos mas adelante
HoltWinters(rain.ts, beta=FALSE, gamma=FALSE)## Holt-Winters exponential smoothing without trend and without seasonal component.
##
## Call:
## HoltWinters(x = rain.ts, beta = FALSE, gamma = FALSE)
##
## Smoothing parameters:
## alpha: 0.02412151
## beta : FALSE
## gamma: FALSE
##
## Coefficients:
## [,1]
## a 24.67819Modelo Holt-Winters con tendencia (doble exponencial)
{r Holt-Winters-tendencia} rm(list=ls(all=TRUE))
use your appropriate directory! setwd(“the directory where you have your data”)
money.data = read.table(“volume-of-money-abs-definition-m.txt”) money.data.ts = ts(money.data[,2],start=c(1960,2) , frequency=12) par(mfrow=c(3,1)) plot(money.data.ts, main=“Time Plot of Volume of Money”) acf(money.data.ts, main=“ACF of Volume of Money”) acf(money.data.ts, type=“partial”, main=“PACF of Volume of Money”)
ajuste manual set up our transformed data and smoothing parameters data <- money.data[,2] N <- length(data) alpha <- 0.7 beta <- 0.5
prepare empty arrays so we can store values forecast <- NULL level <- NULL trend <- NULL
initialize level and trend in a very simple way level[1] <- data [1] trend[1] <- data [2]- data [1]
initialize forecast to get started forecast[1] <- data [1] forecast[2] <- data [2]
loop to build forecasts for( n in 2:N ) { level[n] <- alpha* data [n] +(1-alpha)(level[n-1]+trend[n-1]) trend[n] <- beta(level[n] - level[n-1]) + (1-beta)*trend[n-1] forecast[n+1] <- level[n] + trend[n] }
display your calculated forecast values forecast[3:N]
verificamos con la rutina HoltWinters() m <- HoltWinters(data, alpha = 0.7, beta = 0.5, gamma = FALSE)
plot(m, main=“Holt Winters Fitting of Money Volume with Bogus Parameters”)
usamos la funcion Holt-winters para que busque los alpha y beta optimos m <- HoltWinters(data, gamma = FALSE) plot(m, main=“Holt Winters Fitting of Money Volume with Bogus Parameters”) ```
Modelo Holt-Winters con tendencia y estacionalidad (triple exponencial)
par(mfrow=c(2,1))
plot(AirPassengers,main="Air Passengers time series in thousands")
plot(log10(AirPassengers),main="log10 Air Passengers time series in thousands") #al tomar logaritmos convertimos una estacionalidad multiplicativa por otra aditiva
AirPassengers.SES <- HoltWinters( log10(AirPassengers), beta=FALSE, gamma=FALSE )
AirPassengers.SES$SSE## [1] 0.3065102AirPassengers.SES$alpha #returns 0.9999339## [1] 0.9999339AirPassengers.SES$beta #returns false## [1] FALSEAirPassengers.SES$gamma #returns false## [1] FALSEAirPassengers.HW <- HoltWinters( log10(AirPassengers) ) #con tendencia y estacionalidad
AirPassengers.HW## Holt-Winters exponential smoothing with trend and additive seasonal component.
##
## Call:
## HoltWinters(x = log10(AirPassengers))
##
## Smoothing parameters:
## alpha: 0.326612
## beta : 0.005744246
## gamma: 0.8207255
##
## Coefficients:
## [,1]
## a 2.680598830
## b 0.003900787
## s1 -0.031790733
## s2 -0.061224237
## s3 -0.015941495
## s4 0.006307818
## s5 0.014138008
## s6 0.067260071
## s7 0.127820295
## s8 0.119893006
## s9 0.038321663
## s10 -0.014181699
## s11 -0.085995400
## s12 -0.044672707AirPassengers.HW$SSE## [1] 0.0383026AirPassengers.HW$alpha #returns 0.326612## alpha
## 0.326612AirPassengers.HW$beta #returns 0.005744246## beta
## 0.005744246AirPassengers.HW$gamma #returns 0.8207255## gamma
## 0.8207255#forecast con la funcion Holt-Winters del paquete forecast
rm(list=ls(all=TRUE))
AirPassengers.HW <- HoltWinters(log10(AirPassengers))
AirPassengers.forecast <- forecast(AirPassengers.HW)
AirPassengers.forecast## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Jan 1961 2.652709 2.630898 2.674520 2.619351 2.686066
## Feb 1961 2.627176 2.604218 2.650134 2.592065 2.662287
## Mar 1961 2.676360 2.652297 2.700422 2.639560 2.713160
## Apr 1961 2.702510 2.677380 2.727640 2.664077 2.740942
## May 1961 2.714241 2.688076 2.740406 2.674225 2.754257
## Jun 1961 2.771264 2.744092 2.798436 2.729708 2.812820
## Jul 1961 2.835725 2.807571 2.863878 2.792667 2.878782
## Aug 1961 2.831698 2.802586 2.860811 2.787174 2.876222
## Sep 1961 2.754028 2.723977 2.784079 2.708069 2.799987
## Oct 1961 2.705425 2.674454 2.736396 2.658059 2.752791
## Nov 1961 2.637512 2.605638 2.669386 2.588765 2.686259
## Dec 1961 2.682736 2.649974 2.715497 2.632631 2.732840
## Jan 1962 2.699518 2.661306 2.737731 2.641078 2.757959
## Feb 1962 2.673986 2.635014 2.712957 2.614383 2.733588
## Mar 1962 2.723169 2.683445 2.762894 2.662416 2.783923
## Apr 1962 2.749319 2.708848 2.789790 2.687424 2.811214
## May 1962 2.761050 2.719838 2.802262 2.698022 2.824078
## Jun 1962 2.818073 2.776126 2.860020 2.753921 2.882226
## Jul 1962 2.882534 2.839857 2.925211 2.817265 2.947803
## Aug 1962 2.878508 2.835105 2.921910 2.812129 2.944886
## Sep 1962 2.800837 2.756714 2.844960 2.733356 2.868318
## Oct 1962 2.752234 2.707395 2.797074 2.683658 2.820811
## Nov 1962 2.684322 2.638770 2.729873 2.614656 2.753987
## Dec 1962 2.729545 2.683285 2.775805 2.658796 2.800294par(mfrow=c(1,1))
plot(AirPassengers.forecast, xlim=c(1959, 1963))
AirPassengers.HW## Holt-Winters exponential smoothing with trend and additive seasonal component.
##
## Call:
## HoltWinters(x = log10(AirPassengers))
##
## Smoothing parameters:
## alpha: 0.326612
## beta : 0.005744246
## gamma: 0.8207255
##
## Coefficients:
## [,1]
## a 2.680598830
## b 0.003900787
## s1 -0.031790733
## s2 -0.061224237
## s3 -0.015941495
## s4 0.006307818
## s5 0.014138008
## s6 0.067260071
## s7 0.127820295
## s8 0.119893006
## s9 0.038321663
## s10 -0.014181699
## s11 -0.085995400
## s12 -0.044672707