INTRODUCCIÓN

Se realizará un informe con una base de datos extraidas del BCRP Data https://estadisticas.bcrp.gob.pe/estadisticas/series/mensuales/resultados/PN01803AM/html , que muestra la producción agropecuaria en miles de toneladas de porcinos en un periodo mensual de enero del 2003 hasta noviembre del 2023.

El objetivo del siguiente análisis es:

-Trazar, examinar y preparar series para modelar.

-Extrae el componente de estacionalidad de la serie temporal.

-Probar la estacionalidad y aplicar las transformaciones apropiadas.

-Elija el orden de un modelo ARIMA.

-Pronostica la serie.

PRIMERA ETAPA:

Lectura y preparación de los datos.

PAQUETES

Para iniciar a utlizar los comandos debemos usar las siguientes librerias

library('forecast')
library('tseries')
library('readxl')
#NOTA: si faltase algun comando se debe instalar, usar como ejemplo  install.packages("readxl")
datos=read_excel("C:/Users/Joel/Desktop/Porcino_MT.xlsx")
datosts = ts(datos,start =c(2003,1),frequency=12)
head(datos)
## # A tibble: 6 × 1
##   Porcino_MT
##        <dbl>
## 1       9.94
## 2      10.1 
## 3      10.4 
## 4      10.3 
## 5      10.5 
## 6      10.4
plot.ts(datosts,lwd=3, 
        main ="Serie Orginal de Producción Porcino (Miles de Toneladas)",
        ylab="Miles de Tonaladas",
        xlab="Tiempo ",col="lightblue")

Como podemos observar en la serie, se presenta tendencia y estacionaliedad

par(mfrow=c(1,2))

acf(datosts,main ="ACF de Porcino",)

pacf(datosts,main ="PACF de Porcino")

adf.test(datosts)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  datosts
## Dickey-Fuller = -2.5522, Lag order = 6, p-value = 0.3433
## alternative hypothesis: stationary

Al revisar el test de Dickey Fuller se puede notar una marcada presencia de raiz unitaria debido a que el \(p-valor=0.3433 > 0.05\), por lo que la hipotesis nula del test no se rechaza (\(H_o\) : Presencia de Raiz unitaria )

SEGUNDA ETAPA:

Una vez determinado de nuestra serie original no es apropiada para un análisis estadístico de forma directa se procede a realizar transformaciones para en ese sentido realizar el análisis apropiado, de tal forma que cumplamos todas las condiciones necesarias.

REMOVER LA ESTACIONALIEDAD

Uno de los métodos más utilizados es removiendo la estacionalidad, por ello realizaremos el siguiente análisis.

#install.packages("seasonal")
library(seasonal)
data = seas(datosts)
plot(data, main ="Serie Original VS Serie Ajustada",
     ylab ="Miles de Toneladas" ,
     xlab= "Tiempo"
     )

Se muestra la serie y su grafico ajustado

descomposicion = decompose(datosts)

plot(descomposicion)

Genera un análisis visual de la serie original, tendencia, estacionaliedad y el ruido de la serie

des_por = descomposicion$trend
plot(des_por, main="Tendencia de la Produccion de Porcino(Serie descompuesta)", 
        ylab="Miles de Tonaladas",
        xlab="Tiempo ",col="darkred",lwd=4)

par(mfrow=c(1,2))
acf(na.omit(des_por))
pacf(na.omit(des_por))

adf.test(na.omit(des_por))
## 
##  Augmented Dickey-Fuller Test
## 
## data:  na.omit(des_por)
## Dickey-Fuller = -1.5827, Lag order = 6, p-value = 0.7513
## alternative hypothesis: stationary

Al revisar el test de Dickey Fuller se puede notar una marcada presencia de raiz unitaria debido a que el \(p-valor=0.7513 > 0.05\), por lo que la hipotesis nula del test no se rechaza (\(H_o\) : Presencia de Raiz unitaria )

Para corregir la presencia de raiz unitaria, usaremos un metodo que es la SEGUNDA diferencia

pd_por= diff( des_por)
par(mfrow=c(1,2))
plot(des_por,main = "Serie Original ",
        xlab="Tiempo ",col="darkred",lwd=2)
plot(pd_por, main = "Serie diferenciada ",
        xlab="Tiempo ",col="darkblue",lwd=2)

Comparamos ambas series, y notamos que en la serie diferenciada la media es constante y la varianza semiconstante para verificar si la serie es estacionaria, usariamos el test de Dickey Fuller

par(mfrow=c(1,2))
acf(na.omit(pd_por),main ="ACF - 2da Dif")
pacf(na.omit(pd_por),main ="PACF - 2da Dif")

Asimismo notamos un decaemiento lento y cercano en la franja cercana a cero (\(0\))

adf.test(na.omit(pd_por))
## 
##  Augmented Dickey-Fuller Test
## 
## data:  na.omit(pd_por)
## Dickey-Fuller = -3.8589, Lag order = 6, p-value = 0.01675
## alternative hypothesis: stationary

l revisar el test de Dickey Fuller se puede notar la ausencia de raiz unitaria debido a que el \(p-valor= 0.01675 < 0.05\), por lo que la hipotesis nula del test Se Rechaza (\(H_o\) : Presencia de Raiz unitaria )

TERCERA ETAPA

Ajustamos el modelo.

auto.arima(pd_por, seasonal=FALSE)
## Series: pd_por 
## ARIMA(1,0,1) with non-zero mean 
## 
## Coefficients:
##          ar1     ma1    mean
##       0.7673  0.7547  0.0441
## s.e.  0.0454  0.0771  0.0073
## 
## sigma^2 = 0.0002352:  log likelihood = 656.73
## AIC=-1305.45   AICc=-1305.28   BIC=-1291.56

La representación de una serie ARIMA se expresa en términos de sus componentes autorregresivos (AR), de media móvil (MA), y la diferenciación (I). En este caso, se indica que se tiene un modelo \(ARIMA(0,0,1)\) con media cero, la expresión simplificada sería:

\[Y_t= 0.7673$Y_{t-1}$ + 0.7547\varepsilon_{t-1}+ 0.0441+ \varepsilon_t\]

Esto indica que la serie \(X_t\) en cada periodo \(t\) se compone del término de error blanco en ese periodo y del \(96.41\%\) del término de error blanco en el periodo anterior.

CUARTA ETAPA

Verificación de los parámetros.

Ajustamos el modelo.

modeloarima<-auto.arima(pd_por, seasonal=FALSE)
modeloarima
## Series: pd_por 
## ARIMA(1,0,1) with non-zero mean 
## 
## Coefficients:
##          ar1     ma1    mean
##       0.7673  0.7547  0.0441
## s.e.  0.0454  0.0771  0.0073
## 
## sigma^2 = 0.0002352:  log likelihood = 656.73
## AIC=-1305.45   AICc=-1305.28   BIC=-1291.56
tsdisplay(residuals(modeloarima), lag.max=15, main='(1,0,1) Model Residuals')

Hay un claro patrón presente en ACF que se repite en el retardo 12. Esto sugiere que nuestro modelo puede ser mejor con una especificación diferente como p=12 o q=4.

QUINTA ETAPA

prediccion <- forecast(modeloarima, h=5,level = 95)
prediccion
##          Point Forecast       Lo 95     Hi 95
## Dec 2023     0.04624404 -0.02978676 0.1222748
## Jan 2024     0.04573657 -0.03086505 0.1223382
## Feb 2024     0.04534720 -0.03158847 0.1222829
## Mar 2024     0.04504845 -0.03208320 0.1221801
## Apr 2024     0.04481923 -0.03242756 0.1220660
prediccion <- forecast(modeloarima, h=20,level = 95)
plot(prediccion)