library(writexl)
## Warning: package 'writexl' was built under R version 4.5.1
library(tseries)
## Warning: package 'tseries' was built under R version 4.5.1
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(forecast)
## Warning: package 'forecast' was built under R version 4.5.1
url="https://raw.githubusercontent.com/vmoprojs/DataLectures/master/pib_ec_const.csv"
datos=read.delim(url,sep=";")
pib=datos$PIB
pib_s=ts(pib,start = c(2000,1),frequency = 4)
HoltWinters(pib_s)
## Holt-Winters exponential smoothing with trend and additive seasonal component.
## 
## Call:
## HoltWinters(x = pib_s)
## 
## Smoothing parameters:
##  alpha: 0.886355
##  beta : 0.09704259
##  gamma: 0.105474
## 
## Coefficients:
##           [,1]
## a  16461504.63
## b    -86488.33
## s1   -49526.50
## s2   -38760.76
## s3    16514.29
## s4    46626.38
plot(pib_s)

adf.test(pib_s)
## Warning in adf.test(pib_s): p-value greater than printed p-value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  pib_s
## Dickey-Fuller = 0.22608, Lag order = 4, p-value = 0.99
## alternative hypothesis: stationary
write_xlsx(datos,"Datos_PIB.xlsx")

Al ver el gráfico de la serie, se ve que las fluctuaciones o la estacionalidad son más pronunciadas en los años más recientes, donde la tendencia del PIB es más alta, que en los años iniciales, cuando la tendencia era más baja. Este patrón es un indicador clave de un modelo multiplicativo.

p-valor de 0.99. Dado que el p-valor es mayor a 0.05, se concluye que la serie no es estacionaria.

Los coeficientes estacionales (s1 a s4) indican los valores estacionales de cada trimestre. Los trimestres con los valores más bajos son el trimestre 1 (s1) con -49526.50 y el trimestre 2 (s2) con -38760.76. Los trimestres con los valores más altos son el trimestre 4 (s4) con 46626.38 y el trimestre 3 (s3) con 16514.29.

pib_sl=log(pib_s)
plot(pib_sl)

adf.test(pib_sl)
## Warning in adf.test(pib_sl): p-value greater than printed p-value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  pib_sl
## Dickey-Fuller = 0.75773, Lag order = 4, p-value = 0.99
## alternative hypothesis: stationary

Para estabilizar la varianza, se aplicó logarítmo a la serie, un p-valor de 0.99, lo que confirma que la serie transformada sigue sin ser estacionaria.

pib_sld=diff(pib_sl)
pib_slds=diff(pib_sld,4)
adf.test(pib_slds)
## Warning in adf.test(pib_slds): p-value smaller than printed p-value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  pib_slds
## Dickey-Fuller = -5.4262, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationary

La prueba ADF en la serie diferenciada resultó en p-valor de 0.01 Con un p-valor menor a 0.05, se puede afirmar que la serie es estacionaria.

acf(pib_slds)

pacf(pib_slds)

mod=auto.arima(pib_sl,trace = T)
## 
##  ARIMA(2,2,2)(1,0,1)[4]                    : Inf
##  ARIMA(0,2,0)                              : -373.4734
##  ARIMA(1,2,0)(1,0,0)[4]                    : -400.867
##  ARIMA(0,2,1)(0,0,1)[4]                    : -422.2602
##  ARIMA(0,2,1)                              : -424.4106
##  ARIMA(0,2,1)(1,0,0)[4]                    : -422.26
##  ARIMA(0,2,1)(1,0,1)[4]                    : Inf
##  ARIMA(1,2,1)                              : -422.5483
##  ARIMA(0,2,2)                              : -422.5264
##  ARIMA(1,2,0)                              : -402.8318
##  ARIMA(1,2,2)                              : -420.5152
## 
##  Best model: ARIMA(0,2,1)
mod
## Series: pib_sl 
## ARIMA(0,2,1) 
## 
## Coefficients:
##           ma1
##       -0.9152
## s.e.   0.0452
## 
## sigma^2 = 0.0003178:  log likelihood = 214.28
## AIC=-424.56   AICc=-424.41   BIC=-419.75

La función probó varios modelos y el mejor se identificó como ARIMA(0,2,1) AIC = -424.56, AICC = -424.41 y BIC = -419.75.

\(y''_{t}=-0.9152\epsilon_{t-1}+\epsilon_{t}\)

con un error estandar 0.0452

pib_slp=predict(mod,19)
pib_sp=exp(pib_slp$pred)
pib_sp
##          Qtr1     Qtr2     Qtr3     Qtr4
## 2021 16474318 16416481 16358847 16301415
## 2022 16244185 16187156 16130327 16073697
## 2023 16017267 15961034 15904999 15849161
## 2024 15793518 15738071 15682819 15627760
## 2025 15572895 15518223 15463742

La predicción SARIMA para este trimestre 2025, es de 15,463,742. Para el período total de 2021 a 2025, la predicción muestra una tendencia decreciente, comenzando en 16,474,318 en el primer trimestre de 2021 y finalizando en 15,463,742 en el tercer trimestre de 2025. .

ts.plot(pib_s,pib_sp,col=c("purple","green"))