Resumen

La serie que vamos a utilizar es la del IPC tomada de la pagina del Banco de la Republica, anexa “http://www.banrep.gov.co/es/indicadores-inflacion-basica-y-su-variacion-anual” , con el nombre de IPC, contiene datos en intervalos mensuales durante el mes de Diciembre de 1999 a enero de 2019.

Carga de base de datos:

IPC <- read.csv("IPC.csv")
IPC$mes=paste(IPC$mes,"-01",sep="")
IPC$mes=as.Date(IPC$mes)
IPC <- as.ts(IPC[,2])

A continuacion se presenta el grafico de la serie:

plot(IPC,type="line")

Test de Box-Cox:

Para verificar si la serie necesita una transformacion, utilizaremos este test, tengamos en cuenta que si el test nos da 0 necesita una transformacion logaritmica, si nos da 1 no necesita tranformacion, si es mayor a uno necesita una tranformacion especial.

BoxCox.lambda(IPC,method="guerrero",lower=0)
## [1] 4.102259e-05
BoxCox.lambda(IPC,method="loglik",lower=0)
## [1] 0.25
BoxCox(IPC)

Dado que la prueba Box-Cox al hacer un intervalo de confianza contiene al 0, concluimos que la serie de tiempo requiere una transformacion logaritmica.

#Aplicacion de la transformacion:

lIPC<- log(IPC)
plot(lIPC)

En esta grafica vemos el cambio de la dimension de los datos para que la varianza no se vea tan afectada por la escala de los datos originales.

Revision de ACF Y PACF

acf(lIPC)

pacf(lIPC)

Del ACF y el PACF podemos concluir que la serie tiene un posible AR 1. Ahora vamos a determinar si la serie tiene una ra?z unitaria por lo que ser?a necesario modelarlo como un ARIMA(p,1,q).

adf.test(lIPC)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  lIPC
## Dickey-Fuller = -2.8206, Lag order = 6, p-value = 0.2311
## alternative hypothesis: stationary

La prueba de ra?z unitaria tiene como \(H_{0}\) que la serie contenga alguna raiz unitaria por lo que si no rechazamos dicha hipotesis seria necesario diferenciar).

lIPCdif=diff(lIPC)
plot(lIPCdif)

acf(lIPCdif)

pacf(lIPCdif)

Al diferenciar observamos que su media es igual o cercana a 0 y que muy probablemente su varianza al mantenerse muy parecida en toda la serie, podria ser contante, dado lo anterior podriamos decir que la serie es estacionaria.

Modelo 1

Primero vamos a ajustar un modelo ARMA con la funcion auto.arima

mod1=auto.arima(lIPCdif,max.p =12 ,max.q = 12, start.p = 0, start.q = 0, ic="bic")
mod1
## Series: lIPCdif 
## ARIMA(1,0,0) with zero mean 
## 
## Coefficients:
##          ar1
##       0.3581
## s.e.  0.0628
## 
## sigma^2 estimated as 0.002137:  log likelihood=379.48
## AIC=-754.96   AICc=-754.91   BIC=-748.09

Teniendo en cuenta el modelo anterior y que en el PACF el rezago 12 parece altamente correlacionado vamos a justar un modelo SARIMA, ya que nos indicaria una posible estacionalidad anual en la serie.

lmodeIPC=arima(lIPC,order = c(1,1,0),seasonal = list(order = c(2,1,1), period = 12))
coeftest(lmodeIPC)
## 
## z test of coefficients:
## 
##       Estimate Std. Error  z value  Pr(>|z|)    
## ar1   0.377291   0.064890   5.8143 6.089e-09 ***
## sar1 -0.318313   0.072038  -4.4187 9.931e-06 ***
## sar2 -0.142458   0.070210  -2.0290   0.04246 *  
## sma1 -0.999996   0.055826 -17.9127 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Outliers modelo 1

Ahora vamos a ver cuantos y que tipo de outliers hay en este modelo.

lresi= residuals(lmodeIPC)
plot(lresi)

lcoef= coefs2poly(lmodeIPC)
out=locate.outliers(lresi,lcoef)
out
##   type ind    coefhat     tstat
## 1   AO 176 -0.1024059 -5.237119
## 4   LS 177  0.1406253  4.409050
xreg = outliers.effects(out, 230)
sal1=arima(IPC,order=c(1,1,0),xreg=xreg ,seasonal = list(order = c(2,1,1), period = 12))
coeftest(sal1)
## 
## z test of coefficients:
## 
##        Estimate Std. Error  z value  Pr(>|z|)    
## ar1    0.459370   0.067225   6.8334 8.294e-12 ***
## sar1  -0.454705   0.075837  -5.9958 2.025e-09 ***
## sar2  -0.234448   0.074748  -3.1365   0.00171 ** 
## sma1  -0.999995   0.063397 -15.7736 < 2.2e-16 ***
## AO176 -0.321096   0.204481  -1.5703   0.11635    
## LS177  0.053399   0.340583   0.1568   0.87541    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

En este caso existen dos outliers, uno de cambio de nivel (177) y otro aditivo (176) pero ninguno es significativo para el modelo por lo que no los tendremos en cuenta para el siguiente analisis de residuales.

El modelo final es el que sigue:

sal=arima(IPC,order=c(1,1,0) ,seasonal = list(order = c(2,1,1), period = 12))
coeftest(sal)
## 
## z test of coefficients:
## 
##       Estimate Std. Error  z value  Pr(>|z|)    
## ar1   0.431843   0.067334   6.4134 1.423e-10 ***
## sar1 -0.418486   0.075648  -5.5320 3.166e-08 ***
## sar2 -0.204531   0.074583  -2.7423    0.0061 ** 
## sma1 -0.999995   0.064743 -15.4457 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

##Analisis de residuales para el modelo 1

En esta parte miraremos si el modelo 1 tiene posibles correlaciones.

residuales=sal$residuals
plot(residuales)

acf(residuales)

pacf(residuales)

De lo anterior podemos decir que tiene una autocorrelacion peque?a en el rezago 4, pero no significativa.

A continuacion veremos si los residuales cumplen supuestos de normalidad y de independencia.

jarque.bera.test(residuales)
## 
##  Jarque Bera Test
## 
## data:  residuales
## X-squared = 80.001, df = 2, p-value < 2.2e-16
Box.test(residuales, lag = (length(residuales)/4), type = "Ljung-Box", fitdf = 2)
## 
##  Box-Ljung test
## 
## data:  residuales
## X-squared = 68.5, df = 55.5, p-value = 0.1129

Seg?n el test de Jarque Bera se rechaza la normalidad de los residuales. Y del test de Box-Ljung decimos que a un nivel de significancia del \(5\%\) no se rechaza la independencia.

Dado que se rechaza la normalidad para el test de Jarque Bera observaremos el posible comportamiento distribucional de los residuales.

hist(residuales)

qqnorm(residuales)
qqline(residuales)

Una posible raz?n para que se rechaze la normalidad es que los residuales tienen colas pesadas como se ve en los graficos anteriores ademas de tener una asimetria negativa, podriamos decir que sus posibles distribuciones son una t-student o la distribucion generalizada del error.

Ahora veremos si el modelo podria tener una parte GARCH ya que cabe la posibilidad de que los residuales podrian no ser normales.

Box.test(lmodeIPC$residual^2)
## 
##  Box-Pierce test
## 
## data:  lmodeIPC$residual^2
## X-squared = 2.514, df = 1, p-value = 0.1128

Seg?n el test de Box-Pirce podemos decir que a un nivel de significancia del \(5\%\) se rechaza la hipotesis de que el modelo tenga una parte GARCH.

##Cusum modelo 1

Con los siguentes graficos veremos la estabilidad de los parametros (CUSUM) y la estabilidad de la varianza (CUSUMSQ)

lcum=cumsum(lresi)/sd(lresi)
lN=length(lresi)
lcumq=cumsum(lresi^2)/sum(lresi^2)
Af=0.948 
co=0.14422
lLS=Af*sqrt(lN)+2*Af*c(1:length(lresi))/sqrt(lN)
lLI=-lLS
lLQS=co+(1:length(lresi))/lN
lLQI=-co+(1:length(lresi))/lN
plot(lcum,type="l",ylim=c(min(lLI),max(lLS)),xlab="t",ylab="",main="CUSUM")
lines(lLS,type="S",col="red")
lines(lLI,type="S",col="red")

plot(lcumq,type="l",xlab="t",ylab="",main="CUSUMSQ")        
lines(lLQS,type="S",col="red")                            
lines(lLQI,type="S",col="red")

Podemos ver que la estabilidad de los parametros y de la varianza se tiene en el modelo.

##Pronostico

Ahora veremos los pronosticos de el modelo para los siguientes 12 periodos.

fit= arima(IPC,order=c(1,1,0),include.mean = F,seasonal = list(order = c(2,1,1), period = 12))
pronostico=  forecast(object=fit,h=12)                
plot(pronostico,ylim=c(0,12))

plot(pronostico$residuals,ylim=c(-2,2))