MODELO SARIMA

Al modelo ARIMA se le pueden hacer varias modificaciones para tener en cuenta el comportamiento estacional y no estacionario. A menudo, la dependencia del pasado tiende a ocurrir con más fuerza en múltiplos de algunos rezagos estacionales \(s\)

• Los datos económicos mensuales, hay una fuerte componente anual que ocurre en rezagos que son múltiplos de \(s = 12\), debido a las fuertes conexiones de toda la actividad con el año calendario.

• Es apropiado introducir polinomios de media móvil y autorregresivos. que se identifican con los rezagos estacionales. El autorregresivo estacional puro resultante modelo de promedio móvil, \(ARMA (P, Q)_s\),

\[ \Phi_{P}\left(B^{s}\right) x_{t}=\Theta_{Q}\left(B^{s}\right) w_{t} \]

Donde \[ \Phi_{P}\left(B^{s}\right)=1-\Phi_{1} B^{s}-\Phi_{2} B^{2 s}-\cdots-\Phi_{P} B^{P s} \] y \[ \Theta_{Q}\left(B^{s}\right)=1+\Theta_{1} B^{s}+\Theta_{2} B^{2 s}+\cdots+\Theta_{Q} B^{Q s} \]

son los operador autorregresivos estacional y de media móvil estacional de órdenes \(P\) y \(Q\), respectivamente, con período estacional \(s\).

Ejemplo:

Una Serie de medias móviles autorregresivas estacionales de primer orden podrían escribirse como:

\[\left(1-\Phi B^{12}\right) x_{t}=\left(1+\Theta B^{12}\right) w_{t}\] o \[x_{t}=\Phi x_{t-12}+w_{t}+\Theta w_{t-12}\]

Definición: El modelo autorregresivo estacional multiplicativo integrado de media móvil, o modelo SARIMA, de Box y Jenkins (1970) es dado por

\[\Phi_{P}\left(B^{s}\right) \phi(B) \nabla_{s}^{D} \nabla^{d} x_{t}=\alpha+\Theta_{Q}\left(B^{s}\right) \theta(B) w_{t}\]

donde \(w_t\) es el proceso habitual de ruido blanco. El modelo general se denota como \(ARIMA (p, d, q) × (P, D, Q)_s\).

• Las componentes autorregresivas y de media móvil ordinarias están representados por polinomios \(\phi(B)\) y \(\theta(B)\) de ordenes \(p\) y \(q\), respectivamente.

• Las componentes autorregresivo y de medias moviles estacional \(\Phi(B)\) y \(\Theta(B)\)

• Las diferencias ordinarias y estacionales son respectivamente \(\nabla^{d}=(1-B)^d\) y \(\nabla_{s}^{D}=(1-B^s)^D\)

Ejemplo:

El modelo \(SARIMA(0, 1, 1) × (0, 1, 1)_{12}\) se representa como:

\[\left(1-B^{12}\right)(1-B) x_{t}=\left(1+\Theta B^{12}\right)(1+\theta B) w_{t}\]

Tempdup

Temperatura promedio mensual (en grados Fahrenheit) registrada en Dubuque 1/1964 - 12/1975

Cargar las librerias

library(car)
library(urca)
library(forecast)
library(tseries)
library(FitAR)
library(highcharter)

Cargar los datos

library(readr)
tempdub <- read_csv("tempdub.dat")

Formato serie de tiempo

z1<- ts(tempdub, start = c(1964,1), end=c(1975,12), frequency = 12 )
z1

##       Jan  Feb  Mar  Apr  May  Jun  Jul  Aug  Sep  Oct  Nov  Dec
## 1964 24.7 25.7 30.6 47.5 62.9 68.5 73.7 67.9 61.1 48.5 39.6 20.0
## 1965 16.1 19.1 24.2 45.4 61.3 66.5 72.1 68.4 60.2 50.9 37.4 31.1
## 1966 10.4 21.6 37.4 44.7 53.2 68.0 73.7 68.2 60.7 50.2 37.2 24.6
## 1967 21.5 14.7 35.0 48.3 54.0 68.2 69.6 65.7 60.8 49.1 33.2 26.0
## 1968 19.1 20.6 40.2 50.0 55.3 67.7 70.7 70.3 60.6 50.7 35.8 20.7
## 1969 14.0 24.1 29.4 46.6 58.6 62.2 72.1 71.7 61.9 47.6 34.2 20.4
## 1970  8.4 19.0 31.4 48.7 61.6 68.1 72.2 70.6 62.5 52.7 36.7 23.8
## 1971 11.2 20.0 29.6 47.7 55.8 73.2 68.0 67.1 64.9 57.1 37.6 27.7
## 1972 13.4 17.2 30.8 43.7 62.3 66.4 70.2 71.6 62.1 46.0 32.7 17.3
## 1973 22.5 25.7 42.3 45.2 55.5 68.9 72.3 72.3 62.5 55.6 38.0 20.4
## 1974 17.6 20.5 34.2 49.2 54.8 63.8 74.0 67.1 57.7 50.8 36.8 25.5
## 1975 20.4 19.6 24.6 41.3 61.8 68.5 72.0 71.1 57.3 52.5 40.6 26.2

summary(tempdub)

##     tempdub     
##  Min.   : 8.40  
##  1st Qu.:27.32  
##  Median :48.60  
##  Mean   :46.27  
##  3rd Qu.:63.12  
##  Max.   :74.00

autoplot(z1)#

plot(decompose(z1))

hchart(z1)

Identificación

par(mfrow=c(3,3))
plot(z1,type="o")
acf(z1, lag.max=40)
pacf(z1, lag.max=40)
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), lag.max=40)
pacf(diff(z1), lag.max=40)
plot(diff(z1,12),type="o")
acf(diff(z1,12), lag.max=40)
pacf(diff(z1,12), lag.max=40)

Test DF sobre la serie desestacionalizada

plot(ur.df(diff(z1,12),type="none",lag=12))

summary(ur.df(diff(z1,12),type="none",lag=12))

## 
## ############################################### 
## # 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 
## -11.8377  -2.5662   0.5987   2.5522  10.8391 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## z.lag.1      -1.64348    0.36600  -4.490 1.82e-05 ***
## z.diff.lag1   0.72942    0.33382   2.185  0.03108 *  
## z.diff.lag2   0.63916    0.31374   2.037  0.04412 *  
## z.diff.lag3   0.61574    0.29625   2.078  0.04008 *  
## z.diff.lag4   0.50539    0.27730   1.823  0.07119 .  
## z.diff.lag5   0.46527    0.25846   1.800  0.07468 .  
## z.diff.lag6   0.56632    0.23824   2.377  0.01924 *  
## z.diff.lag7   0.52325    0.21942   2.385  0.01887 *  
## z.diff.lag8   0.49019    0.19540   2.509  0.01363 *  
## z.diff.lag9   0.44101    0.17154   2.571  0.01153 *  
## z.diff.lag10  0.45063    0.14953   3.014  0.00323 ** 
## z.diff.lag11  0.32356    0.12630   2.562  0.01182 *  
## z.diff.lag12 -0.11132    0.09405  -1.184  0.23925    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.263 on 106 degrees of freedom
## Multiple R-squared:  0.5952, Adjusted R-squared:  0.5456 
## F-statistic: 11.99 on 13 and 106 DF,  p-value: 1.379e-15
## 
## 
## Value of test-statistic is: -4.4904 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau1 -2.58 -1.95 -1.62

Dickey Fuller del paquete tseries

adf.test(diff(z1,12))

## Warning in adf.test(diff(z1, 12)): p-value smaller than printed p-value

## 
##  Augmented Dickey-Fuller Test
## 
## data:  diff(z1, 12)
## Dickey-Fuller = -4.9185, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationary

test phillips- perron

plot(ur.pp(diff(z1,12),type="Z-tau",
     model="constant", lags="long"))

summary(ur.pp(diff(z1,12),type="Z-tau",
     model="constant", lags="long"))

## 
## ################################## 
## # Phillips-Perron Unit Root Test # 
## ################################## 
## 
## Test regression with intercept 
## 
## 
## Call:
## lm(formula = y ~ y.l1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -11.0779  -3.1268   0.3738   2.7058  13.0120 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.03669    0.41033  -0.089    0.929
## y.l1         0.08989    0.08660   1.038    0.301
## 
## Residual standard error: 4.695 on 129 degrees of freedom
## Multiple R-squared:  0.008283,   Adjusted R-squared:  0.000595 
## F-statistic: 1.077 on 1 and 129 DF,  p-value: 0.3012
## 
## 
## Value of test-statistic, type: Z-tau  is: -10.579 
## 
##          aux. Z statistics
## Z-tau-mu           -0.0858
## 
## Critical values for Z statistics: 
##                      1pct      5pct     10pct
## critical values -3.480998 -2.883488 -2.578338

Como la serie desestacionalizada es estacionaria vamos a plantear algunos modelos

par(mfrow=c(1,3))
plot(diff(z1,12),type="o")
acf(diff(z1,12), lag.max=40)
pacf(diff(z1,12), lag.max=40)

Ajuste del Modelo

modelo1<-stats::arima(z1,
        order=c(0,0,11), 
        seasonal=list(order=c(0,1,1),
        period=12), fixed = c(0,0,0,0,0,NA,0,0,0,0,NA,NA)) 
modelo1

## 
## Call:
## stats::arima(x = z1, order = c(0, 0, 11), seasonal = list(order = c(0, 1, 1), 
##     period = 12), fixed = c(0, 0, 0, 0, 0, NA, 0, 0, 0, 0, NA, NA))
## 
## Coefficients:
##       ma1  ma2  ma3  ma4  ma5     ma6  ma7  ma8  ma9  ma10     ma11     sma1
##         0    0    0    0    0  0.1560    0    0    0     0  -0.1521  -1.0000
## s.e.    0    0    0    0    0  0.0904    0    0    0     0   0.0946   0.1033
## 
## sigma^2 estimated as 11.17:  log likelihood = -361.81,  aic = 731.63

tt <- modelo1$coef[which(modelo1$coef!=0)]/sqrt(diag(modelo1$var.coef))
1 - pt(abs(tt),(modelo1$nobs - length(modelo1$coef[which(modelo1$coef!=0)])))

##        ma6       ma11       sma1 
## 0.04339467 0.05511757 0.00000000

BIC(modelo1)

## [1] 743.1567

Diagnóstico

et<-residuals(modelo1)
x1.fit <- fitted(modelo1)

par(mfrow=c(3,2))
plot(z1,type="l",lty=2)
lines(x1.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] 111  27

acf(abs(et)) #Mide Estructura Heteroced?stica

tsdiag(modelo1)

Test de Autocorrelacion de Ljung-Box

\(H_0: r1=r2=r3=...=rlag=0\)

\(H_a\): Al menos una es dif de cero

Box.test(et,lag=30,type="Ljung-Box")

## 
##  Box-Ljung test
## 
## data:  et
## X-squared = 25.288, df = 30, p-value = 0.7109

Test de Normalidad basado en Sesgo y Curtosis

\(H_0:\) Datos provienen de una Dist. Normal

\(H_a:\) Los datos no provienen de una Dist. Normal

jarque.bera.test(et)

## 
##  Jarque Bera Test
## 
## data:  et
## X-squared = 1.5898, df = 2, p-value = 0.4516

Test de Aleatoriedad

\(H_0:\) Residuales exhiben un comport. de Aleatoriedad \(H_a:\) Residuales no exhiben estructura (Tendencia, o cualquier otro comportamiento predecible)

runs.test(as.factor(sign(et)), alternative="two.sided")

## 
##  Runs Test
## 
## data:  as.factor(sign(et))
## Standard Normal = -2.0051, p-value = 0.04495
## alternative hypothesis: two.sided

Pronóstico Fuera Muestra

plot(forecast(modelo1,h=10, fan=T))
lines(fitted(modelo1), col="red")