Modelos de media móvil autorregresiva (ARMA)
Al combinar los dos modelos MA y AR, obtenemos lo que se llama un modelo de promedio móvil autoregresivo (ARMA).
El caso más simple, es el proceso ARMA (\(1\),\(1\)) como
\[\begin{equation} y_{t}=\phi y_{t-1}+\varepsilon_{t}+\theta \varepsilon_{t-1} \tag{3.1} \end{equation}\]
En la discusión relacionada con los procesos de promedio móvil, notamos que el ACF daría una indicación del orden del proceso, mientras que el PACF disminuiría lentamente. En contraste, cuando discutimos los modelos autoregresivos, notamos que el PACF daría una indicación del orden del proceso, mientras que el ACF disminuiría lentamente.
Cuando consideramos el proceso ARMA combinado, generalmente notamos que ambas funciones deberían decaer ligeramente y, como tal, puede ser difícil descifrar el orden del modelo ARMA combinado.
Un ejemplo de las funciones ACF y PACF para un modelo ARMA (\(1\),\(1\)). A pesar del hecho de que sabemos que estamos tratando con un modelo ARMA (1,1), las autocorrelaciones en ACF y PACF parecen diferir de cero en
Algunas librerias
library(foreign)
library(ggfortify)## Warning: package 'ggfortify' was built under R version 4.1.2## Warning: package 'ggplot2' was built under R version 4.1.2library(forecast)## Warning: package 'forecast' was built under R version 4.1.2library(gridExtra)## Warning: package 'gridExtra' was built under R version 4.1.2library(seasonal)## Warning: package 'seasonal' was built under R version 4.1.2library(lattice)
library(zoo)## Warning: package 'zoo' was built under R version 4.1.2library(urca)## Warning: package 'urca' was built under R version 4.1.2library(dynlm)## Warning: package 'dynlm' was built under R version 4.1.2Simulación ARMA(\(1,1\))
set.seed(123)
arma1 = arima.sim(n = 1000, list(ar = 0.5,ma=0.7))
ggtsdisplay(arma1)
set.seed(123)
arma1 = arima.sim(n = 500, list(ar = 0.5,ma=0.5))
ggtsdisplay(arma1)
Los modelos ARMA (incluidos los términos AR y MA) tienen ACF y PACF que ambos se reducen a 0. Estos son los más complicados porque el orden no será particularmente obvio. Básicamente, solo tiene que adivinar que uno o dos términos de cada tipo pueden ser necesarios y luego ver qué sucede cuando estima el modelo.
ARMA(p,q)
\[\begin{equation} y_{t}=\phi_1 y_{t-1}+ \phi_2 y_{t-2}+ \cdots\phi_p y_{t-p} + \varepsilon_{t}+\theta_1 \varepsilon_{t-1}+ \theta_2 \varepsilon_{t-2}+ \cdots+ \theta_q \varepsilon_{t-q} \end{equation}\]
Diferenciación
A menudo, la diferenciación se usa para dar cuenta de la no estacionariedad que ocurre en forma de tendencia y / o estacionalidad.
La diferencia \(x_t-x_{t-1}\) puede expresarse como \((1-B)x_{t}\).
Una notación alternativa para una diferencia es
\[(\boldsymbol{1-B})=\triangle\]
Así \[\triangle x_t = \boldsymbol{(1-B)}x_t = x_t-x_{t-1}\]
Un superindice define una diferencia de retraso igual al subíndice. Por ejemplo,
\[\triangle^{12} x_t = (\boldsymbol{1-B}^{12})x_t = x_t-x_{t-12}\]
Este tipo de diferencia a menudo se usa con datos mensuales que exhiben estacionalidad
data("AirPassengers")
autoplot(AirPassengers)
AirPassengers## Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 1949 112 118 132 129 121 135 148 148 136 119 104 118
## 1950 115 126 141 135 125 149 170 170 158 133 114 140
## 1951 145 150 178 163 172 178 199 199 184 162 146 166
## 1952 171 180 193 181 183 218 230 242 209 191 172 194
## 1953 196 196 236 235 229 243 264 272 237 211 180 201
## 1954 204 188 235 227 234 264 302 293 259 229 203 229
## 1955 242 233 267 269 270 315 364 347 312 274 237 278
## 1956 284 277 317 313 318 374 413 405 355 306 271 306
## 1957 315 301 356 348 355 422 465 467 404 347 305 336
## 1958 340 318 362 348 363 435 491 505 404 359 310 337
## 1959 360 342 406 396 420 472 548 559 463 407 362 405
## 1960 417 391 419 461 472 535 622 606 508 461 390 432diff(AirPassengers)## Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 1949 6 14 -3 -8 14 13 0 -12 -17 -15 14
## 1950 -3 11 15 -6 -10 24 21 0 -12 -25 -19 26
## 1951 5 5 28 -15 9 6 21 0 -15 -22 -16 20
## 1952 5 9 13 -12 2 35 12 12 -33 -18 -19 22
## 1953 2 0 40 -1 -6 14 21 8 -35 -26 -31 21
## 1954 3 -16 47 -8 7 30 38 -9 -34 -30 -26 26
## 1955 13 -9 34 2 1 45 49 -17 -35 -38 -37 41
## 1956 6 -7 40 -4 5 56 39 -8 -50 -49 -35 35
## 1957 9 -14 55 -8 7 67 43 2 -63 -57 -42 31
## 1958 4 -22 44 -14 15 72 56 14 -101 -45 -49 27
## 1959 23 -18 64 -10 24 52 76 11 -96 -56 -45 43
## 1960 12 -26 28 42 11 63 87 -16 -98 -47 -71 42autoplot(diff(AirPassengers))
ggtsdisplay(AirPassengers)
ggtsdisplay(diff(AirPassengers))
ggtsdisplay(diff(AirPassengers,12))
En general, se sabe que el punto de partida para el análisis en series de tiempo estacionarias. Es decir, que los datos involucrados cumplan unas exigencias mínimas:
- Ser constante en media (No tener tendencia en ninguna forma)
- Ser constante en varianza.
Si bien el criterio visual es importante para determinar el tratamiento a tomar, son necesarios criterios técnicos para tomar tales decisiones. Es así que se postulan las
De una forma muy general y sin entrar en detalles, notaremos como \(I(0)\) a las Series integradas de orden 0. Dichas series son aquellas que no presentan problemas de estacionalidad. Análogamente, notaremos como \(I(1)\) a las series no estacionarías.
Así, en el caso que sea de nuestro interés saber cuándo una serie es bien comportada en el sentido de estacionariedad, tenemos la prueba aumentada de Dickey-Fuller que establecerá si la serie es integrada de orden 1. Es decir, formalmente \[
\begin{cases}
H_0:y_t\sim I(1)\\
H_1:y_t\sim I(0)
\end{cases}
\]
En la práctica, se utiliza la prueba de Dickey-Fuller aumentada con los comandos en ,
- \(\texttt{ur.df()}\) y
- \(\texttt{adf.test()}\)
de los paquetes \(\texttt{urca}\) y \(\texttt{tseries}\), respectivamente. Ahora, mientras que la prueba \(\texttt{adf.test()}\) reporta directamente el p-valor de la prueba de raíz unitaria, la función \(\texttt{ut.test()}\) reporta los valores de la estadística de prueba y de los valores críticos con los valores de significancia \(1\%, 5\%\) y \(10\%\).
FASES
La metodología de Box y Jenkins se resume en cuatro fases:La primera fase consiste en identificar el posible modelo ARIMA que sigue la serie, lo que requiere:
- Decidir qué transformaciones aplicar para convertir la serie observada en una serie estacionaria.
- Determinar un modelo ARMA para la serie estacionaria, es decir, los órdenes p y q de su estructura autorregresiva y de media móvil.
La segunda fase: estimación, En esta etapa se estiman los coeficientes de los términos autorregresivos y de media móvil incluidos en el modelo, cuyo número de rezagos \(p\) y \(q\) ya han sido identificados en la etapa anterior.
La tercera fase es el diagnostico, donde se comprueba que los residuos no tienen estructura de dependencia y siguen un proceso de ruido blanco.
La cuarta fase es la predicción, una vez que se ha obtenido un modelo adecuado se realizan predicciones con el mismo
Identificación de un modelo
Decidir si \(x_t\) necesita ser transformada para eliminar la no estacionariedad en media o la no estacionariedad en varianza (heteroscedasticidad). Puede ser conveniente utilizar logaritmos de la serie o aplicar la transformación de Box‐Cox.
Determinación del grado \(d\) de diferenciación adecuado.
En general, la falta de estacionariedad se manifiesta en que los coeficientes de la función de autocorrelación estimada tienden a decrecer muy lentamente.
- Decidir los valores de \((p, q)\), y si existe una componente estacional, decidir los órdenes de los operadores estacionales \((P, Q)\). Para este apartado se utilizan las funciones de autocorrelación (ACF) y autocorrelación parcial (PACF) según el siguiente cuadro:
Metodología Box Jenkins
Caffeine
library(car)## Warning: package 'car' was built under R version 4.1.2## Warning: package 'carData' was built under R version 4.1.2library(urca)
library(forecast)
library(tseries)## Warning: package 'tseries' was built under R version 4.1.2library(FitAR)## Warning: package 'FitAR' was built under R version 4.1.2## Warning: package 'leaps' was built under R version 4.1.2## Warning: package 'bestglm' was built under R version 4.1.2Cargamos los datos de Caffeine almacenados en R
data(Caffeine) \(X\) concebtracion de la cafeina en un proceso industrial
- Encabezado, seis primeros datos
head(Caffeine)## Time Series:
## Start = 1
## End = 6
## Frequency = 1
## [1] 0.429 0.443 0.451 0.455 0.440 0.433- Longitud de la serie de tiempo
length(Caffeine) ## [1] 178z1<-Caffeine # escribir en una nueva variable
head(z1)## Time Series:
## Start = 1
## End = 6
## Frequency = 1
## [1] 0.429 0.443 0.451 0.455 0.440 0.433summary(z1)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.3590 0.3930 0.4040 0.4055 0.4170 0.4550Identificación
par(mfrow=c(1,3))
plot(z1)
acf(z1)
pacf(z1)
par(mfrow=c(2,3)) # garaficar en 2 filas y tres columnas
plot(z1,type="o")
acf(z1)
pacf(z1)
plot(diff(z1),type="o")
abline(h=2*sqrt(var(diff(z1))),col="red",lty=2)
abline(h=-2*sqrt(var(diff(z1))),col="red",lty=2)
acf(diff(z1))
pacf(diff(z1))
Dickey-Fuller aumentado
adf.test(z1)##
## Augmented Dickey-Fuller Test
##
## data: z1
## Dickey-Fuller = -2.2021, Lag order = 5, p-value = 0.4918
## alternative hypothesis: stationaryPaquete urca
plot(ur.df(z1,type="none",lag=10))
summary(ur.df(z1,type="none",lag=10))##
## ###############################################
## # Augmented Dickey-Fuller Test Unit Root Test #
## ###############################################
##
## Test regression none
##
##
## Call:
## lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.016334 -0.003355 -0.000041 0.003139 0.016869
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## z.lag.1 -0.001127 0.001092 -1.032 0.303836
## z.diff.lag1 0.259361 0.074479 3.482 0.000645 ***
## z.diff.lag2 -0.153953 0.076086 -2.023 0.044737 *
## z.diff.lag3 -0.143746 0.077017 -1.866 0.063860 .
## z.diff.lag4 0.190497 0.077316 2.464 0.014830 *
## z.diff.lag5 -0.752542 0.076461 -9.842 < 2e-16 ***
## z.diff.lag6 0.282672 0.076556 3.692 0.000307 ***
## z.diff.lag7 0.017966 0.078367 0.229 0.818971
## z.diff.lag8 -0.069246 0.077084 -0.898 0.370397
## z.diff.lag9 0.169848 0.076747 2.213 0.028344 *
## z.diff.lag10 -0.363719 0.074531 -4.880 2.6e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.005649 on 156 degrees of freedom
## Multiple R-squared: 0.4677, Adjusted R-squared: 0.4302
## F-statistic: 12.46 on 11 and 156 DF, p-value: < 2.2e-16
##
##
## Value of test-statistic is: -1.0316
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau1 -2.58 -1.95 -1.62Phillips Perron
pp.test(z1)## Warning in pp.test(z1): p-value smaller than printed p-value##
## Phillips-Perron Unit Root Test
##
## data: z1
## Dickey-Fuller Z(alpha) = -30.755, Truncation lag parameter = 4, p-value
## = 0.01
## alternative hypothesis: stationaryPruebas de estacionariedad para la serie diferenciada
adf.test(diff(z1))## Warning in adf.test(diff(z1)): p-value smaller than printed p-value##
## Augmented Dickey-Fuller Test
##
## data: diff(z1)
## Dickey-Fuller = -7.0074, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationaryp-value = 0.01 La serie es estacionaria
pp.test(diff(z1))## Warning in pp.test(diff(z1)): p-value smaller than printed p-value##
## Phillips-Perron Unit Root Test
##
## data: diff(z1)
## Dickey-Fuller Z(alpha) = -120.99, Truncation lag parameter = 4, p-value
## = 0.01
## alternative hypothesis: stationaryAjuste del Modelo
modelo1<-stats::arima(z1,
order=c(0,1,5), fixed=c(NA,NA,0,0,NA))
modelo1##
## Call:
## stats::arima(x = z1, order = c(0, 1, 5), fixed = c(NA, NA, 0, 0, NA))
##
## Coefficients:
## ma1 ma2 ma3 ma4 ma5
## 0.1550 0.0543 0 0 -0.7928
## s.e. 0.0616 0.0432 0 0 0.0550
##
## sigma^2 estimated as 3.225e-05: log likelihood = 661.32, aic = -1314.63tt <- modelo1$coef[which(modelo1$coef!=0)]/sqrt(diag(modelo1$var.coef))
tt## ma1 ma2 ma5
## 2.514177 1.255769 -14.4228281 - pt(abs(tt),(modelo1$nobs - length(modelo1$coef[which(modelo1$coef!=0)])))## ma1 ma2 ma5
## 0.006419305 0.105441679 0.000000000BIC(modelo1)## [1] -1301.927Diagnóstico
Calcula y almacena los residuales
et<-residuals(modelo1)
et## Time Series:
## Start = 1
## End = 178
## Frequency = 1
## [1] 4.289996e-04 1.088110e-02 5.169201e-03 2.267882e-03 -1.184523e-02
## [6] -3.531189e-03 -2.199546e-03 -6.628547e-03 7.623142e-04 6.885014e-03
## [11] 6.012180e-03 2.218805e-03 -1.597031e-04 -2.152766e-03 -5.748008e-04
## [16] -6.061131e-03 3.208374e-03 3.547250e-03 2.641616e-03 -4.920557e-03
## [21] -2.786956e-03 -5.720643e-03 -3.059108e-04 -1.708727e-03 -1.421995e-02
## [26] 2.667552e-03 1.934567e-03 2.301033e-03 -1.866908e-02 -4.386787e-03
## [31] -4.186510e-03 -6.321693e-04 -2.107038e-03 5.669918e-03 -4.119073e-03
## [36] 9.911246e-04 -1.032524e-02 -1.129313e-03 -7.807010e-03 9.879667e-04
## [41] -7.977644e-03 -9.855343e-04 -3.344547e-03 3.415883e-03 5.358028e-03
## [46] 4.717696e-03 -2.778695e-03 5.267381e-04 7.883740e-04 -3.871638e-03
## [51] 1.223569e-03 2.831824e-03 9.302720e-04 3.291706e-03 -3.585331e-03
## [56] 3.339924e-04 3.821478e-04 6.744869e-04 -2.531383e-03 4.502059e-03
## [61] -6.259247e-03 -3.983884e-03 4.681208e-04 -2.825816e-03 -5.282848e-05
## [66] -5.785574e-03 1.733222e-03 -3.576173e-03 -5.769675e-03 3.013915e-03
## [71] 5.278836e-03 -2.600711e-03 -6.726098e-03 -4.386180e-03 4.374231e-04
## [76] 1.133641e-02 6.163495e-03 2.100117e-03 6.854865e-03 7.186071e-03
## [81] -2.491893e-03 -1.129779e-03 1.696340e-02 8.909602e-04 -7.368773e-03
## [86] 1.101269e-03 -2.639948e-03 -2.144464e-04 -1.247546e-04 9.190384e-03
## [91] 4.604946e-04 -4.661159e-03 1.519860e-03 1.926882e-03 3.901122e-03
## [96] 3.652248e-03 -3.466753e-03 7.537083e-03 1.553073e-03 1.443560e-03
## [101] -4.415807e-03 -1.139323e-03 -3.609316e-03 -5.146559e-03 2.133028e-03
## [106] 2.450627e-03 -5.396867e-03 1.837376e-03 -2.068067e-03 -2.087446e-03
## [111] -1.624586e-03 -9.112043e-04 -4.314041e-03 -4.921789e-03 -2.658447e-03
## [116] -4.607638e-03 -8.654597e-04 4.963349e-03 3.377412e-03 -3.899942e-03
## [121] -1.623046e-02 -3.961078e-03 4.429787e-03 -7.927234e-04 2.787652e-03
## [126] -9.253249e-03 1.141761e-03 -5.162279e-03 -6.889658e-03 -1.744305e-02
## [131] -2.259105e-03 1.620081e-02 1.752096e-02 -2.056207e-03 -3.461007e-03
## [136] 2.857728e-03 -1.410223e-03 -9.047805e-03 8.846828e-03 1.337745e-02
## [141] -1.286716e-03 2.353410e-03 -1.046583e-02 -3.492810e-03 -4.286140e-03
## [146] -1.165850e-03 -2.721585e-03 5.188169e-03 2.575096e-03 -2.078280e-03
## [151] 3.257479e-03 4.510739e-04 8.660942e-04 -6.117297e-03 7.253025e-03
## [156] -2.090717e-04 9.961800e-04 2.543169e-03 -8.296932e-03 4.896777e-03
## [161] 5.262118e-04 1.442433e-03 -9.236137e-03 -4.224458e-03 -4.962024e-03
## [166] 6.415346e-03 1.041862e-02 -1.284357e-03 -1.071569e-02 7.964071e-04
## [171] -1.455547e-03 -7.558337e-03 2.316795e-04 9.879718e-03 8.806599e-05
## [176] -6.704064e-03 -2.957907e-03 -4.993782e-03fitted(modelo1)## Time Series:
## Start = 1
## End = 178
## Frequency = 1
## [1] 0.4285710 0.4321189 0.4458308 0.4527321 0.4518452 0.4365312 0.4251995
## [8] 0.4186285 0.4102377 0.4191150 0.4299878 0.4387812 0.4461597 0.4451528
## [15] 0.4375748 0.4320611 0.4237916 0.4274527 0.4333584 0.4369206 0.4357870
## [22] 0.4297206 0.4203059 0.4177087 0.4192199 0.4053324 0.4120654 0.4146990
## [29] 0.4186691 0.4083868 0.4001865 0.3936322 0.3911070 0.4033301 0.4131191
## [36] 0.4120089 0.4133252 0.4031293 0.3968070 0.3910120 0.3909776 0.3899855
## [43] 0.3893445 0.3915841 0.3946420 0.4072823 0.4137787 0.4134733 0.4112116
## [50] 0.4078716 0.3997764 0.4031682 0.4060697 0.4067083 0.4135853 0.4086660
## [57] 0.4066179 0.4063255 0.4045314 0.4044979 0.4092592 0.4019839 0.3965319
## [64] 0.3988258 0.3920528 0.3967856 0.3932668 0.3945762 0.3927697 0.3859861
## [71] 0.3937212 0.3986007 0.3987261 0.3953862 0.3875626 0.3836636 0.3988365
## [78] 0.4118999 0.4181451 0.4258139 0.4254919 0.4181298 0.4150366 0.4291090
## [85] 0.4253688 0.4188987 0.4206399 0.4042144 0.4031248 0.4088096 0.4185395
## [92] 0.4216612 0.4164801 0.4180731 0.4130989 0.4173478 0.4254668 0.4204629
## [99] 0.4274469 0.4265564 0.4254158 0.4231393 0.4156093 0.4101466 0.4028670
## [106] 0.4085494 0.4123969 0.4091626 0.4150681 0.4110874 0.4066246 0.4089112
## [113] 0.4063140 0.4029218 0.3986584 0.3966076 0.3918655 0.3940367 0.4036226
## [120] 0.4098999 0.4092305 0.3909611 0.3815702 0.3837927 0.3862123 0.4022532
## [127] 0.3948582 0.3921623 0.3868897 0.3764431 0.3632591 0.3587992 0.3814790
## [134] 0.4080562 0.4204610 0.4181423 0.4084102 0.3930478 0.3841532 0.3966226
## [141] 0.4102867 0.4106466 0.4204658 0.4014928 0.3862861 0.3821658 0.3787216
## [148] 0.3838118 0.3924249 0.3990783 0.3977425 0.4035489 0.4001339 0.3991173
## [155] 0.3937470 0.3992091 0.3990038 0.3994568 0.4072969 0.3921032 0.3974738
## [162] 0.3975576 0.3972361 0.3932245 0.3839620 0.3775847 0.3835814 0.4032844
## [169] 0.4057157 0.3972036 0.3924555 0.3825583 0.3747683 0.3831203 0.3939119
## [176] 0.3957041 0.3939579 0.3899938z1.fit<-z1-etpar(mfrow=c(3,2))
plot(z1,type="l",lty=2)
lines(z1.fit,type="l",col="red")
plot(scale(et),type="l",main="Residuales")
abline(h=2*sqrt(var(scale(et))),col="red",lty=2)
abline(h=-2*sqrt(var(scale(et))),col="red",lty=2)
acf(et)
pacf(et)
qqPlot(scale(et))## [1] 29 133acf(abs(et)) #Mide Estructura Heterocedástica
tsdiag(modelo1) #Resumen Gráfico (Residuales Estandarizados, ACF, Box Test) 
Test de Autocorrelacion de Ljung-Box
\[\begin{align} H_o:& r1=r2=r3=...=rlag=0\\ H_a&: \text{Al menos una es dif de cero} \end{align}\]
\(H_0:\) Los datos se distribuyen de forma independiente (es decir, las correlaciones en la población de la que se toma la muestra son 0, de modo que cualquier correlación observada en los datos es el resultado de la aleatoriedad del proceso de muestreo).
\(H_a:\) Los datos no se distribuyen de forma independiente.
OBJETIVO: “NO RECHAZAR \(H_0\)”
Box.test(et,lag=30,type="Ljung-Box")##
## Box-Ljung test
##
## data: et
## X-squared = 42.6, df = 30, p-value = 0.06352Test de Normalidad basado en Sesgo y Curtosis
\[\begin{align} H_0&: \text{Datos provienen de una Dist. Normal}\\ H_a&: \text{Los datos no provienen de una Dist. Normal} \end{align}\]
jarque.bera.test(et)##
## Jarque Bera Test
##
## data: et
## X-squared = 17.778, df = 2, p-value = 0.0001379p-value = 0.0001379 Los residuales no son normales
Test de Aleatoriedad
\[\begin{align} H_0:& \text{Residuales exhiben un comport. de Aleatoriedad}\\ H_a:& \text{Residuales no exhiben estructura (Tendencia, o cualquier otro comportamiento predecible)} \end{align}\]
OBJETIVO “NO RECHAZAR \(H_0\)
runs.test(as.factor(sign(et)), alternative="two.sided")##
## Runs Test
##
## data: as.factor(sign(et))
## Standard Normal = -1.7512, p-value = 0.07992
## alternative hypothesis: two.sidedp-value = 0.07992 no se rechaza. Los residuales exben un comportamiento de aleatoriedad
Pronóstico Fuera Muestra
plot(forecast(modelo1,h=10, fan=T))
lines(fitted(modelo1), col="red")
#legend(legend=c("Original","Ajustado"),col=c("black","red"))
pred4<-forecast(modelo1,h=10)
pred4## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 179 0.3762336 0.3689558 0.3835113 0.3651031 0.3873640
## 180 0.3758927 0.3647743 0.3870110 0.3588887 0.3928966
## 181 0.3812073 0.3670275 0.3953871 0.3595211 0.4028934
## 182 0.3835522 0.3668633 0.4002410 0.3580288 0.4090756
## 183 0.3875110 0.3686439 0.4063781 0.3586562 0.4163657
## 184 0.3875110 0.3684019 0.4066200 0.3582862 0.4167357
## 185 0.3875110 0.3681630 0.4068589 0.3579208 0.4171011
## 186 0.3875110 0.3679270 0.4070949 0.3575599 0.4174620
## 187 0.3875110 0.3676938 0.4073281 0.3572033 0.4178187
## 188 0.3875110 0.3674634 0.4075586 0.3568508 0.4181712Ejemplo 2
library(forecast)
library(tseries)
library(car)
library(urca)
library(FitAR)1. SERIESB2
Precios de las acciones de IBM, Precio de cierre de las acciones ordinarias, diario, 29 de junio de 1959 al 30 de junio de 1960
data(SeriesB2)
plot(SeriesB2)
summary(SeriesB2)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 401.0 416.5 428.0 438.3 450.0 541.0boxplot(SeriesB2)
Como serie de tiempo
SeriesB2 <- as.ts(SeriesB2)par(mfrow=c(2,3))
plot(SeriesB2)
acf(SeriesB2)
pacf(SeriesB2)
plot(diff(SeriesB2))
acf(diff(SeriesB2))
pacf(diff(SeriesB2))
plot(ur.df(SeriesB2, lag=9))
summary(ur.df(SeriesB2, lag=9))##
## ###############################################
## # Augmented Dickey-Fuller Test Unit Root Test #
## ###############################################
##
## Test regression none
##
##
## Call:
## lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)
##
## Residuals:
## Min 1Q Median 3Q Max
## -13.0017 -2.9336 -0.1673 2.2738 23.4645
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## z.lag.1 0.0007396 0.0007282 1.016 0.3109
## z.diff.lag1 -0.0561974 0.0627986 -0.895 0.3718
## z.diff.lag2 -0.0378272 0.0626570 -0.604 0.5466
## z.diff.lag3 0.0461067 0.0626593 0.736 0.4626
## z.diff.lag4 0.0252789 0.0622310 0.406 0.6850
## z.diff.lag5 -0.0318005 0.0628649 -0.506 0.6134
## z.diff.lag6 0.1314211 0.0631873 2.080 0.0386 *
## z.diff.lag7 0.0385439 0.0637408 0.605 0.5460
## z.diff.lag8 0.0959809 0.0640020 1.500 0.1350
## z.diff.lag9 -0.1323317 0.0642606 -2.059 0.0406 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.906 on 235 degrees of freedom
## Multiple R-squared: 0.0595, Adjusted R-squared: 0.01948
## F-statistic: 1.487 on 10 and 235 DF, p-value: 0.145
##
##
## Value of test-statistic is: 1.0156
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau1 -2.58 -1.95 -1.62plot(ur.df(diff(SeriesB2), lag=9))
summary(ur.df(diff(SeriesB2), lag=9))##
## ###############################################
## # Augmented Dickey-Fuller Test Unit Root Test #
## ###############################################
##
## Test regression none
##
##
## Call:
## lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)
##
## Residuals:
## Min 1Q Median 3Q Max
## -12.6429 -2.6574 0.1096 2.3942 23.7932
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## z.lag.1 -0.85876 0.19702 -4.359 1.96e-05 ***
## z.diff.lag1 -0.17472 0.18856 -0.927 0.355
## z.diff.lag2 -0.21349 0.18188 -1.174 0.242
## z.diff.lag3 -0.16368 0.17344 -0.944 0.346
## z.diff.lag4 -0.13861 0.16332 -0.849 0.397
## z.diff.lag5 -0.16767 0.14925 -1.123 0.262
## z.diff.lag6 -0.03575 0.13579 -0.263 0.793
## z.diff.lag7 0.00513 0.11866 0.043 0.966
## z.diff.lag8 0.10408 0.09522 1.093 0.275
## z.diff.lag9 -0.02613 0.06495 -0.402 0.688
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.914 on 234 degrees of freedom
## Multiple R-squared: 0.5503, Adjusted R-squared: 0.5311
## F-statistic: 28.64 on 10 and 234 DF, p-value: < 2.2e-16
##
##
## Value of test-statistic is: -4.3588
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau1 -2.58 -1.95 -1.62Ajuste del Modelo
modelo1 <- stats::arima((SeriesB2),order=c(0,1,6),
fixed=c(NA,NA,NA,NA,NA,NA))
modelo1##
## Call:
## stats::arima(x = (SeriesB2), order = c(0, 1, 6), fixed = c(NA, NA, NA, NA, NA,
## NA))
##
## Coefficients:
## ma1 ma2 ma3 ma4 ma5 ma6
## -0.0758 -0.0424 0.0808 -0.0029 -0.0684 0.1861
## s.e. 0.0627 0.0637 0.0650 0.0592 0.0611 0.0676
##
## sigma^2 estimated as 24.39: log likelihood = -766.19, aic = 1546.37BIC(modelo1)## [1] 1571.135tt <- modelo1$coef[which(modelo1$coef!=0)] / sqrt(diag(modelo1$var.coef))
1 - pt(abs(tt), (modelo1$nobs - 1))## ma1 ma2 ma3 ma4 ma5 ma6
## 0.113988948 0.253215222 0.107622562 0.480323271 0.132227027 0.003172317Diagnóstico
et <- residuals(modelo1)
yfit <- fitted(modelo1)
par(mfrow=c(3,2))
plot(yfit, type="l", col="red", lty=2, lwd=2)
lines((SeriesB2), col="blue")
legend(10,550,legend=c("Fit","SeriesB2"), col=c("red","blue"), lty=c(2,1), lwd=c(2,1), bty="n")
plot(et);abline(h=mean(et),col="blue",lty=2);abline(h=-2*sd(et), col="red", lty=2);abline(h=2*sd(et), col="red", lty=2)
qqPlot(et)## [1] 229 10acf(et)
pacf(et)
acf(abs(et)) # mide heterocedasticidad
Pruebas al modelo
Box.test(et,lag=24,type="Ljung-Box")##
## Box-Ljung test
##
## data: et
## X-squared = 11.049, df = 24, p-value = 0.9887jarque.bera.test(et)##
## Jarque Bera Test
##
## data: et
## X-squared = 104.27, df = 2, p-value < 2.2e-16runs.test(as.factor(sign(et)))##
## Runs Test
##
## data: as.factor(sign(et))
## Standard Normal = -0.18216, p-value = 0.8555
## alternative hypothesis: two.sidedSi los residuales ressponden a un generador de datos diferente a la dist Normal
require(MASS)## Loading required package: MASSfitdistr(et, "t")## Warning in dt((x - m)/s, df, log = TRUE): NaNs produced## Warning in log(s): Se han producido NaNs## Warning in dt((x - m)/s, df, log = TRUE): NaNs produced## m s df
## 0.02396954 3.37433714 3.39937397
## (0.25524607) (0.28888467) (0.84175942)qqPlot(et, dist="t", df=3.37)
## [1] 229 10Pronóstico
plot(forecast(modelo1,h=4, fan=T), main="Pronóstico SerieB2")
lines((SeriesB2), col="red", lty=2)
library(highcharter)## Warning: package 'highcharter' was built under R version 4.1.2Para visualizar mejor
hchart(forecast(modelo1,h=4, fan=T))