La media móvil integrada autorregresiva estacional, SARIMA o ARIMA estacional, es una extensión de ARIMA que admite explícitamente datos de series temporales univariadas con un componente estacional.
library(knitr)
library(MASS)
library(parallel)
library(mlogit)
library(dplyr)
library(forecast)
library(tseries)
library(highcharter)
Utilizaermos los datos que descargamos del siguiente enlace: https://estadisticas.bcrp.gob.pe/estadisticas/series/trimestrales/resultados/PN02971FQ/html En este caso solo eh traido los datos sin fecha y lo eh guardado localemente como Trimestral9.csv, posteriormente creare las fechas, no olvidar que la data esta en fechas trimestrales.
#Cargamos los datos
desembolso <- read.csv2("Trimestral9.csv")
Los datos son de tipo data.frame debemos convertirlos en series temporales ts y le asignaremos fecha para poder trabajar.start la fecha de inicio y el mes frecuency = 4 en este caso son datos trimestrales
#transformamos los datos en una serie temporal
desembolso_ts <- ts(desembolso, start = c(1999,1),frequency = 4)
desembolso_ts
## Qtr1 Qtr2 Qtr3 Qtr4
## 1999 171.80600 417.90200 210.63400 436.73500
## 2000 553.24700 223.45700 388.84800 319.19300
## 2001 225.63500 503.50300 338.09500 276.36400
## 2002 1553.98600 272.10200 393.44300 682.81100
## 2003 859.65600 81.17600 251.44000 968.82500
## 2004 269.43600 658.02500 251.96800 1355.23100
## 2005 635.39000 74.75800 940.89900 1004.98500
## 2006 53.52590 53.27525 66.00778 436.30462
## 2007 2340.92884 56.79291 305.03845 680.81590
## 2008 203.60268 183.44755 247.43092 531.76492
## 2009 1188.58962 89.05037 1213.08401 738.09087
## 2010 232.78189 1847.08931 243.02061 1938.19287
## 2011 303.41664 101.95346 217.56347 366.85922
## 2012 979.38044 56.74539 109.42919 303.88212
## 2013 559.02118 263.58323 317.23219 137.53266
## 2014 99.74862 1244.03727 702.35217 876.12804
## 2015 943.94783 41.24877 2763.62978 1441.21724
## 2016 1266.29062 140.18375 149.09438 552.48636
## 2017 671.92737 2083.10982 123.10366 166.97152
## 2018 107.01067 80.63800 62.05709 1550.32772
## 2019 519.21696 815.90288 50.22926 477.59026
## 2020 78.49915 3323.14540 2443.07055 4132.67598
## 2021 6120.80278 452.48077 2167.57468 5220.20924
## 2022 119.86665
El test de DICKEY-FULLER (ADF) es una prueba de raíz unitaria para una muestra de una serie de tiempo,nos permite saber si existe estacionalidad.
#Test de DICKEY-FULLER
adf.test(desembolso_ts)
##
## Augmented Dickey-Fuller Test
##
## data: desembolso_ts
## Dickey-Fuller = -2.6458, Lag order = 4, p-value = 0.3099
## alternative hypothesis: stationary
Segun el test tenemos un P-value = 0.3090, lo que indica que no es suficiente para recahzar la hipotesis nula.En conclusion no estacionalidad.
#comande para saber cuantas diferenciacaiones se requieren para un SARIMA
nsdiffs(desembolso_ts)
## [1] 0
La función ACF es usada para identificar el proceso de media móvil (MA) en un modelo ARIMA; mientras que la función PACF se usa para identificar los valores de la parte del proceso autoregresivo (AR).
#funcion de autocorrelacion (ACF)
plot(acf(desembolso_ts))
#funcion de autocorrelacion Parcial (PACF)
plot(pacf(desembolso_ts))
En R tenemos una funcion auto.arima, el cual nos ofrece una alternaiva no unica para nuestro modelo.
#la funcion auto.arima() calcula el mejor modelo ARIMA(p,d,q) de acuerdo a direfentes de criterios AIC,AICC o bic VALUE
model <-auto.arima(desembolso_ts, stepwise = FALSE, approximation = FALSE)
model
## Series: desembolso_ts
## ARIMA(2,1,0)(1,0,0)[4]
##
## Coefficients:
## ar1 ar2 sar1
## -0.7378 -0.4888 -0.2942
## s.e. 0.0960 0.1019 0.1129
##
## sigma^2 = 890746: log likelihood = -759.89
## AIC=1527.78 AICc=1528.24 BIC=1537.86
En conjunto, la notación para un modelo SARIMA se especifica como: SARIMA(p,d,q)(P,D,Q)m — ARIMA(2,1,0)(1,0,0)[4]
checkresiduals(model)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(2,1,0)(1,0,0)[4]
## Q* = 2.7828, df = 5, p-value = 0.7334
##
## Model df: 3. Total lags used: 8
El test nos indica si nuestro modelo esta bien ajustado y si se cumple los requisitos de ruido blanco.
En nuestro modelo1 segun el test tenemos un p-value= 0.7 mayor que el nivel de significancia entonces se concluye que es ruido blanco.
En nuetra primera grafica observamos presencia de ruido blanco, teniendo una media
que gira al rededor del cero.
Observamos que en nuestro modelo1 no hay autocorrelacion, esto quiere decir que un valor determinado no depende de valores anteriores.
adf.test(resid(model))
## Warning in adf.test(resid(model)): p-value smaller than printed p-value
##
## Augmented Dickey-Fuller Test
##
## data: resid(model)
## Dickey-Fuller = -4.8335, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationary
Observamos en nuestro test un p-value = 0.01 menor que nuestro nivel de significancia
confirmamos ruido blanco.
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 2022 Q2 3383.729 2174.2084 4593.249 1533.9272 5233.530
## 2022 Q3 3908.757 2658.3455 5159.169 1996.4175 5821.097
## 2022 Q4 1470.743 162.6141 2778.871 -529.8673 3471.353
## 2023 Q1 3604.116 2085.6974 5122.534 1281.8952 5926.337
## 2023 Q2 2929.729 1408.2145 4451.243 602.7734 5256.684
model2<-auto.arima(desembolso_ts)
model2
## Series: desembolso_ts
## ARIMA(0,1,3)(1,0,0)[4]
##
## Coefficients:
## ma1 ma2 ma3 sar1
## -0.8111 0.0737 0.2243 -0.4858
## s.e. 0.1086 0.1388 0.1169 0.1355
##
## sigma^2 = 887357: log likelihood = -759.29
## AIC=1528.57 AICc=1529.27 BIC=1541.18
checkresiduals(model2)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(0,1,3)(1,0,0)[4]
## Q* = 1.1018, df = 4, p-value = 0.894
##
## Model df: 4. Total lags used: 8
predic2<- forecast(model,h=5)
predic2
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 2022 Q2 3383.729 2174.2084 4593.249 1533.9272 5233.530
## 2022 Q3 3908.757 2658.3455 5159.169 1996.4175 5821.097
## 2022 Q4 1470.743 162.6141 2778.871 -529.8673 3471.353
## 2023 Q1 3604.116 2085.6974 5122.534 1281.8952 5926.337
## 2023 Q2 2929.729 1408.2145 4451.243 602.7734 5256.684
hchart(predic2)