Paquetes

library(TSA)
## 
## Attaching package: 'TSA'
## The following objects are masked from 'package:stats':
## 
##     acf, arima
## The following object is masked from 'package:utils':
## 
##     tar
library(ggfortify)
## Loading required package: ggplot2
library(forecast)
library(fpp)
## Loading required package: fma
## Loading required package: expsmooth
## Loading required package: lmtest
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## Loading required package: tseries
library(seasonal)
library(fpp2)
## 
## Attaching package: 'fpp2'
## The following objects are masked from 'package:fpp':
## 
##     ausair, ausbeer, austa, austourists, debitcards, departures,
##     elecequip, euretail, guinearice, oil, sunspotarea, usmelec
library(gridExtra)
library(urca)
library(foreign)
library(vars)
## Loading required package: MASS
## 
## Attaching package: 'MASS'
## The following objects are masked from 'package:fma':
## 
##     cement, housing, petrol
## Loading required package: strucchange
## Loading required package: sandwich
library(zoo)
library(dynlm)
library(lmtest)
library(strucchange)
library(readxl)

1. Breve descripcion de la serie elegida, asi como su frecuencia, unidad de medida, fuente, etc.

Serie: PM

2006-2018, incluye serie original, desestacionalizada.

Unidad de medida:indice.

Periodicidad: Mensual

Fuente:(www.inegi.org.mx)

2. Graficar los datos, analizar patrones y observaciones atipicas. Apoye su analisis con una descomposicion clasica.

base<-read_xls("C:/Users/rodar/OneDrive/Escritorio/SERIES/INDUSTRIA MANUFACTURERA.xls")
PM<-ts(base$PM, frequency = 12, start=c(2005,1))
autoplot(PM) + xlab("Year")+ylab("INDICE") +
  ggtitle('INDUSTRIA MANUFACTURERA (PM ORIGINAL)')

fit <- decompose(PM,type = "additive")
autoplot(fit)

3. Si es necesario, utilizar transformacion Box-Cox/logaritmos para estabilizar la varianza.

PMlog <-log(PM)
plot(PMlog)

Se utilizo la aplicacion de un logaritmo ya que se hiba a estabilizar la varianza para visualizar el comportamiento de la serie.

4. En caso de estacionalidad, aplicar diferencias estacionales.

PMdif <-diff(PMlog)
plot(PMlog)

PMdif2 <-diff(log(PM))
BoxCox.ar(PM)
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1

## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1

## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1

## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1

## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1

## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1

## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1

## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1

## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1

## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1

## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1

## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1

## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1

## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1

plot(PMdif2)

Sabemos que contiene un valor =0 muy cerca de su centro por lo que sugiere este una transformacion logaritmica.

plot(PM^1.5, type="o")

plot(diff(PM^1.5), type="o")

plot(diff(PMdif2))

Se agrearon dos diferencias.

5. Usar preba Dickey-Fuller para evaluar el orden de integracion de la serie. Diferenciar hasta que la serie sea estacionaria.

summary(ur.df(PM))
## 
## ############################################### 
## # 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 
## -8.6093 -1.1080 -0.2844  2.0807  4.9716 
## 
## Coefficients:
##            Estimate Std. Error t value Pr(>|t|)
## z.lag.1    0.003696   0.003588   1.030    0.308
## z.diff.lag 0.115448   0.138306   0.835    0.408
## 
## Residual standard error: 2.675 on 49 degrees of freedom
## Multiple R-squared:  0.04397,    Adjusted R-squared:  0.004953 
## F-statistic: 1.127 on 2 and 49 DF,  p-value: 0.3323
## 
## 
## Value of test-statistic is: 1.03 
## 
## Critical values for test statistics: 
##      1pct  5pct 10pct
## tau1 -2.6 -1.95 -1.61

No es estacionaria debido a que el valor de P es de 0.3323, y no es apto para 0.5

summary(ur.df(PMdif2))
## 
## ############################################### 
## # 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 
## -0.088854 -0.008099  0.000222  0.022948  0.052844 
## 
## Coefficients:
##            Estimate Std. Error t value Pr(>|t|)    
## z.lag.1    -0.84361    0.18071  -4.668 2.47e-05 ***
## z.diff.lag  0.04486    0.13834   0.324    0.747    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.02737 on 48 degrees of freedom
## Multiple R-squared:  0.405,  Adjusted R-squared:  0.3802 
## F-statistic: 16.33 on 2 and 48 DF,  p-value: 3.884e-06
## 
## 
## Value of test-statistic is: -4.6684 
## 
## Critical values for test statistics: 
##      1pct  5pct 10pct
## tau1 -2.6 -1.95 -1.61

Como podemos ver ya es estacionaria.

6. Usar herramientas ACF, PACF, EACF y criterios de Akaike/ Bayes para construir propuestas de modelos.

ggAcf(PMdif2)

ggPacf(PMdif2) 

eacf(PMdif2)
## AR/MA
##   0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 o o o o o o o o o o o  o  o  o 
## 1 x o o o o o o o o o o  o  o  o 
## 2 x x o o o o o o o o o  o  o  o 
## 3 x o x o o o o o o o o  o  o  o 
## 4 x o x o o o o o o o o  o  o  o 
## 5 o o o o o o o o o o o  o  o  o 
## 6 x o o x o o o o o o o  o  o  o 
## 7 x o o x o o o o o o o  o  o  o
ggtsdisplay(PMdif2)

PM2 <-armasubsets(y=PMdif2,nar= 12, nma=12, y.name="test", ar.method="ols")
## Warning in leaps.setup(x, y, wt = wt, nbest = nbest, nvmax = nvmax,
## force.in = force.in, : 2 linear dependencies found
plot(PM2)

propuesta1<- Arima(PMdif2, order = c(1,0,0),method = "ML")
propuesta1
## Series: PMdif2 
## ARIMA(1,0,0) with non-zero mean 
## 
## Coefficients:
##          ar1    mean
##       0.1540  0.0056
## s.e.  0.1402  0.0044
## 
## sigma^2 estimated as 0.0007603:  log likelihood=113.95
## AIC=-221.9   AICc=-221.4   BIC=-216.04
propuesta2 <-Arima(PMdif2, order = c(0,0,2),method = "ML")
propuesta2
## Series: PMdif2 
## ARIMA(0,0,2) with non-zero mean 
## 
## Coefficients:
##          ma1     ma2    mean
##       0.1863  0.0439  0.0057
## s.e.  0.1784  0.1484  0.0046
## 
## sigma^2 estimated as 0.0007736:  log likelihood=114.02
## AIC=-220.03   AICc=-219.18   BIC=-212.23

7. Como parte del diagnostico del modelo, analizar los residuos y aplicar la prueba Ljung-Box.

tsdiag(propuesta1,gof=12, omit.initial=F)

tsdiag(propuesta2,gof=12, omit.initial=F)

Como podemos ver. en las dos se ajustan los residuos, ya que los valores son arriba de P.

8. Usar la funcion auto.arima () y compare sus resultados con el inciso anterior.

propuesta.auto<-auto.arima(PMdif2)
propuesta.auto
## Series: PMdif2 
## ARIMA(0,0,0) with zero mean 
## 
## sigma^2 estimated as 0.0007783:  log likelihood=112.33
## AIC=-222.67   AICc=-222.59   BIC=-220.72
checkresiduals(propuesta.auto)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,0,0) with zero mean
## Q* = 11.041, df = 10.4, p-value = 0.3883
## 
## Model df: 0.   Total lags used: 10.4
propuesta2 <- auto.arima(PMdif2, stepwise = FALSE, approximation = FALSE)
propuesta2
## Series: PMdif2 
## ARIMA(1,0,2) with zero mean 
## 
## Coefficients:
##           ar1     ma1     ma2
##       -0.8872  1.2367  0.4829
## s.e.   0.3140  0.1993  0.2908
## 
## sigma^2 estimated as 0.0006855:  log likelihood=116.81
## AIC=-225.62   AICc=-224.77   BIC=-217.82
checkresiduals(propuesta2)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(1,0,2) with zero mean
## Q* = 9.4314, df = 7.4, p-value = 0.2556
## 
## Model df: 3.   Total lags used: 10.4

Se puede observar que en las dos propuestas,la media ya es estable asi como los residuales.

9. Presente la ecuacion final y describa.

\[ Yt= -0.8872e_1+1.2367e_{t-1}+0.4829e_{t-2}+et \]

10. Crear un pronostico a dos años.

pronostico<- plot(forecast(propuesta.auto, h=24))

pronostico
## $mean
##      Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 2009                       0   0   0   0   0   0   0
## 2010   0   0   0   0   0   0   0   0   0   0   0   0
## 2011   0   0   0   0   0                            
## 
## $lower
##                  80%         95%
## Jun 2009 -0.03575187 -0.05467777
## Jul 2009 -0.03575187 -0.05467777
## Aug 2009 -0.03575187 -0.05467777
## Sep 2009 -0.03575187 -0.05467777
## Oct 2009 -0.03575187 -0.05467777
## Nov 2009 -0.03575187 -0.05467777
## Dec 2009 -0.03575187 -0.05467777
## Jan 2010 -0.03575187 -0.05467777
## Feb 2010 -0.03575187 -0.05467777
## Mar 2010 -0.03575187 -0.05467777
## Apr 2010 -0.03575187 -0.05467777
## May 2010 -0.03575187 -0.05467777
## Jun 2010 -0.03575187 -0.05467777
## Jul 2010 -0.03575187 -0.05467777
## Aug 2010 -0.03575187 -0.05467777
## Sep 2010 -0.03575187 -0.05467777
## Oct 2010 -0.03575187 -0.05467777
## Nov 2010 -0.03575187 -0.05467777
## Dec 2010 -0.03575187 -0.05467777
## Jan 2011 -0.03575187 -0.05467777
## Feb 2011 -0.03575187 -0.05467777
## Mar 2011 -0.03575187 -0.05467777
## Apr 2011 -0.03575187 -0.05467777
## May 2011 -0.03575187 -0.05467777
## 
## $upper
##                 80%        95%
## Jun 2009 0.03575187 0.05467777
## Jul 2009 0.03575187 0.05467777
## Aug 2009 0.03575187 0.05467777
## Sep 2009 0.03575187 0.05467777
## Oct 2009 0.03575187 0.05467777
## Nov 2009 0.03575187 0.05467777
## Dec 2009 0.03575187 0.05467777
## Jan 2010 0.03575187 0.05467777
## Feb 2010 0.03575187 0.05467777
## Mar 2010 0.03575187 0.05467777
## Apr 2010 0.03575187 0.05467777
## May 2010 0.03575187 0.05467777
## Jun 2010 0.03575187 0.05467777
## Jul 2010 0.03575187 0.05467777
## Aug 2010 0.03575187 0.05467777
## Sep 2010 0.03575187 0.05467777
## Oct 2010 0.03575187 0.05467777
## Nov 2010 0.03575187 0.05467777
## Dec 2010 0.03575187 0.05467777
## Jan 2011 0.03575187 0.05467777
## Feb 2011 0.03575187 0.05467777
## Mar 2011 0.03575187 0.05467777
## Apr 2011 0.03575187 0.05467777
## May 2011 0.03575187 0.05467777