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)
library(rmarkdown)

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

Serie: PG

2006-2018, incluye serie original, desestacionalizada.

Unidad de medida: Millones de pies cubicos.

Periodicidad: Mensual

Fuente: PEMEX. Dirección Corporativa de Finanzas. (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/PRODUCCION DE GAS.xls")
PGN <-ts(base$PG, frequency = 12, start=c(2006,1))
autoplot(PGN) + xlab("Year")+ylab("Millones de pies cubicos") +
  ggtitle('PRODUCCION NACIONAL DE GAS (PGN ORIGINAL)')

#Graficando para encontrar estacionalidad#
ggseasonplot(PGN)

ggseasonplot(PGN, polar = TRUE)

#Nos apoyamos de una descomposicion clasica#
fit<-decompose(PGN)
autoplot(fit)

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

#No aplica la transformacion Box-Cox#

logPGN<- log(PGN)
plot(logPGN)

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

4. En caso de estacionalidad, aplicar diferencias estacionales.

PGNdif<-diff(logPGN)
plot(PGNdif)

Demuestra estacionalidad la serie por lo que se aplico la diferencia para que la media se estabilizara.

BoxCox.ar(PGN)

## 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

## 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

## 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

## 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

## 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

## 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

## 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

## 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

## 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

## 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

## 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

## 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

## 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

## 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

PGNdiflog<-diff(logPGN)

plot(PGNdiflog)

?? ?? sabemos que contiene un valor= 0 muy cerca de su centro por lo que sugiere este una transformacion logaritmica. el maximo esta al rededor 1.5

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

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

plot(diff(PGNdiflog))

Se agregaron dos diferencias con logaritmo para estabilizar la serie, aqui se puede observar que se elevo al 1.5 para tener mejoras en la varianza como en la media.

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

summary(ur.df(PGN))

## 
## ############################################### 
## # 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 
## -30279  -3129   1791   5072  13343 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## z.lag.1    -0.0003578  0.0034483  -0.104    0.917    
## z.diff.lag -0.5865141  0.0673642  -8.707 7.52e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7949 on 141 degrees of freedom
## Multiple R-squared:  0.3499, Adjusted R-squared:  0.3407 
## F-statistic: 37.94 on 2 and 141 DF,  p-value: 6.535e-14
## 
## 
## Value of test-statistic is: -0.1038 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau1 -2.58 -1.95 -1.62

Como podemos ver, nuestro valor de P= 6.535e-14 que nos indica un valor mayor a 0.5

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

ggAcf(PGNdiflog)

ggPacf(PGNdiflog)

eacf(PGNdiflog)

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

ggtsdisplay(PGNdiflog)

con respecto a la pruebas de correlacion y correlacion parcial, podemos observar que se pueden aplicar los siguientes modelos que serian un ar(2) o en su defecto un ma(2).

CRITERIO DE BAYES

PGN.2 <-armasubsets(y=PGNdiflog,nar= 12, nma=12, y.name="test", ar.method="ols")

plot(PGN.2)

se puede aplicar una ar(2) o un ma(2)

propuesta <- Arima(PGNdiflog, order = c(2,1,2),method = "ML")
propuesta

## Series: PGNdiflog 
## ARIMA(2,1,2) 
## 
## Coefficients:
##           ar1     ar2      ma1     ma2
##       -0.0538  0.1592  -1.7584  0.7932
## s.e.   0.1422  0.1294   0.1080  0.1020
## 
## sigma^2 estimated as 0.0017:  log likelihood=252.57
## AIC=-495.14   AICc=-494.7   BIC=-480.32

En la aplicacion de un ar(2) podemos ver que en el criterio de akaike es -497.7 y esto no esta diciendo que tan bien estan ajustados los datos para este modelo de ar.

propuesta2 <- Arima(PGNdiflog, order = c(1,1,1),method = "ML")
propuesta2

## Series: PGNdiflog 
## ARIMA(1,1,1) 
## 
## Coefficients:
##           ar1      ma1
##       -0.5993  -0.9672
## s.e.   0.0682   0.0208
## 
## sigma^2 estimated as 0.001933:  log likelihood=242.83
## AIC=-479.66   AICc=-479.48   BIC=-470.77

De las dos propuestas se escoge el AR(2) por que la prueba de akakike ajustada es menor que el de la propuesta de AR(1) entonces se aprecia que la prueba de akaike ajustada tiene un valor de -479.48 y nos esta diciendo que en relacion con el AR(2) los datos se estan ajustando de una mejor manera ya que es mas pequeño el valor en el AR(1).

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

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

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

Se ajustan mas los residuos en la propuesta dos que en la uno, porque en la pruebas de la ljung-box y la ACF se ven la variaciones mas estables, los valores de p son arriba de 0.5.

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

propuesta.auto<-auto.arima(PGNdiflog)
propuesta.auto

## Series: PGNdiflog 
## ARIMA(1,0,1)(0,1,1)[12] with drift 
## 
## Coefficients:
##          ar1      ma1     sma1   drift
##       0.3541  -0.6215  -0.8771  -1e-04
## s.e.  0.1867   0.1500   0.1474   1e-04
## 
## sigma^2 estimated as 0.0006228:  log likelihood=293.3
## AIC=-576.59   AICc=-576.12   BIC=-562.18

checkresiduals(propuesta.auto)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(1,0,1)(0,1,1)[12] with drift
## Q* = 14.063, df = 20, p-value = 0.8273
## 
## Model df: 4.   Total lags used: 24

propuest.auto2 <- auto.arima(PGNdiflog, stepwise = FALSE, approximation = FALSE)
propuest.auto2

## Series: PGNdiflog 
## ARIMA(1,0,1)(0,1,1)[12] with drift 
## 
## Coefficients:
##          ar1      ma1     sma1   drift
##       0.3541  -0.6215  -0.8771  -1e-04
## s.e.  0.1867   0.1500   0.1474   1e-04
## 
## sigma^2 estimated as 0.0006228:  log likelihood=293.3
## AIC=-576.59   AICc=-576.12   BIC=-562.18

checkresiduals(propuest.auto2)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(1,0,1)(0,1,1)[12] with drift
## Q* = 14.063, df = 20, p-value = 0.8273
## 
## Model df: 4.   Total lags used: 24

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.3541 e_1-0.6215e_{t-1}+et \] Tiene correlacion.

10. Crear un pronostico a dos años.

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

pronostico

## $mean
##              Jan         Feb         Mar         Apr         May
## 2018             -0.11331309  0.07759195 -0.05026112  0.01834441
## 2019 -0.01045285 -0.10360933  0.07994103 -0.05051619  0.01716726
## 2020 -0.01213535                                                
##              Jun         Jul         Aug         Sep         Oct
## 2018 -0.03545813  0.02092041 -0.01227701 -0.06848575  0.03743088
## 2019 -0.03696178  0.01930116 -0.01393720 -0.07016043  0.03575107
## 2020                                                            
##              Nov         Dec
## 2018 -0.04045739  0.02432303
## 2019 -0.04213901  0.02264076
## 2020                        
## 
## $lower
##                   80%         95%
## Feb 2018 -0.145463820 -0.16248339
## Mar 2018  0.044314856  0.02669903
## Apr 2018 -0.083676747 -0.10136591
## May 2018 -0.015088543 -0.03278688
## Jun 2018 -0.068893260 -0.08659275
## Jul 2018 -0.012514994 -0.03021463
## Aug 2018 -0.045712448 -0.06341210
## Sep 2018 -0.101921188 -0.11962084
## Oct 2018  0.003995450 -0.01370420
## Nov 2018 -0.073892777 -0.09159240
## Dec 2018 -0.009112022 -0.02681147
## Jan 2019 -0.043885209 -0.06158323
## Feb 2019 -0.137271828 -0.15509168
## Mar 2019  0.046259503  0.02842958
## Apr 2019 -0.084200099 -0.10203128
## May 2019 -0.016516947 -0.03434829
## Jun 2019 -0.070646021 -0.08847738
## Jul 2019 -0.014383092 -0.03221446
## Aug 2019 -0.047621445 -0.06545281
## Sep 2019 -0.103844675 -0.12167604
## Oct 2019  0.002066832 -0.01576453
## Nov 2019 -0.075823213 -0.09365455
## Dec 2019 -0.011043103 -0.02887426
## Jan 2020 -0.045816538 -0.06364628
## 
## $upper
##                   80%           95%
## Feb 2018 -0.081162368 -0.0641428013
## Mar 2018  0.110869034  0.1284848615
## Apr 2018 -0.016845493  0.0008436717
## May 2018  0.051777368  0.0694757053
## Jun 2018 -0.002023005  0.0156764822
## Jul 2018  0.054355805  0.0720554364
## Aug 2018  0.021158420  0.0388580691
## Sep 2018 -0.035050314 -0.0173506624
## Oct 2018  0.070866315  0.0885659639
## Nov 2018 -0.007021997  0.0106776294
## Dec 2018  0.057758082  0.0754575296
## Jan 2019  0.022979500  0.0406775197
## Feb 2019 -0.069946838 -0.0521269890
## Mar 2019  0.113622548  0.1314524688
## Apr 2019 -0.016832285  0.0009988988
## May 2019  0.050851466  0.0686828072
## Jun 2019 -0.003277533  0.0145538286
## Jul 2019  0.052985405  0.0708167693
## Aug 2019  0.019747053  0.0375784173
## Sep 2019 -0.036476178 -0.0186448143
## Oct 2019  0.069435318  0.0872666793
## Nov 2019 -0.008454811  0.0093765281
## Dec 2019  0.056324627  0.0741557887
## Jan 2020  0.021545838  0.0393755814

11. Utilizando métodos de pronósticos simples y/o suavizamiento exponencial generar un pronóstico para dos años. Elegir el método que presente la raíz del error medio al cuadrado más pequeño.

data(PGN)

## Warning in data(PGN): data set 'PGN' not found

autoplot(PGN)+
 ggtitle("PRODUCCION NACIONAL DE GAS")+
 ylab("PRODUCCION")+xlab("AÑO")

Datos establecidos a partir de 2006-2018.

PGN3<-window(PGN,start = 2006,end = c(2017,12))

Trazando algunos pronosticos.

autoplot(PGN3)+
  forecast::autolayer(meanf(PGN3, h=11), PI=FALSE, series="Mean")+
  forecast::autolayer(naive(PGN3, h=11), PI=FALSE, series="Naïve")+
  forecast::autolayer(snaive(PGN3, h=11), PI=FALSE, series= "Seasonal naïve")+
  ggtitle("Pronosticos para la produccion de gas nacional")+
  xlab("Año")+ylab("Produccion")+
  guides(colour= guide_legend(title="Pronostico"))

Utilizamos metodo de pronostico simple.

PGNfit1<- meanf(PGN3,h=10)
PGNfit2<- rwf(PGN3,h=10)
PGNfit3<- snaive(PGN3,h=10)
autoplot(window(PGN, start=2006))+
  forecast::autolayer(PGNfit1, series="Mean", PI=FALSE)+
  forecast::autolayer(PGNfit2, series="Naïve", PI=FALSE)+
  forecast::autolayer(PGNfit3, series= "Seasonal naïve", PI=FALSE)+
  ggtitle("Pronosticos para la produccion de gas nacional")+
  xlab("Año")+ylab("Produccion")+
  guides(colour= guide_legend(title="Pronostico"))

PGN4<- window(PGN, start=2018)
accuracy(PGNfit1, PGN4)

##                         ME     RMSE      MAE        MPE      MAPE     MASE
## Training set  4.547474e-12 19788.11 15674.74  -1.180777  8.607752 1.289316
## Test set     -3.935890e+04 39358.90 39358.90 -25.857952 25.857952 3.237440
##                  ACF1
## Training set 0.849851
## Test set           NA

accuracy(PGNfit2, PGN4)

##                      ME     RMSE      MAE        MPE     MAPE      MASE
## Training set  -54.38874 9870.156 8045.443 -0.1802689 4.288304 0.6617725
## Test set     3074.01487 3074.015 3074.015  2.0195616 2.019562 0.2528510
##                    ACF1
## Training set -0.5877821
## Test set             NA

accuracy(PGNfit3, PGN4)

##                       ME     RMSE      MAE       MPE     MAPE     MASE
## Training set   -796.6549 15396.11 12157.42 -0.869997 6.553287 1.000000
## Test set     -12879.4706 12879.47 12879.47 -8.461535 8.461535 1.059392
##                   ACF1
## Training set 0.9222148
## Test set            NA

12. Analizar un VAR, necesita incluir los puntos siguientes.

rm(list = ls()) # limpiar la consola

base2<-read_xls("C:/Users/rodar/OneDrive/Escritorio/SERIES/PNG Y LIQ.xls")

Declaramos ambas variables como series de tiempo.

ts.P1 <- ts(base2$P, start = c(2006,1), frequency = 12)
ts.Q2 <- ts(base2$Q, start = c(2006,1), frequency = 12)

autoplot(ts.P1)

autoplot(ts.Q2)

ggseasonplot(ts.P1)

ggseasonplot(ts.Q2)

mes. <- season(ts.P1)
reg.P <- lm(ts.P1 ~ time(ts.P1) + mes.)
summary(reg.P)

## 
## Call:
## lm(formula = ts.P1 ~ time(ts.P1) + mes.)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -48999 -11060   5760  13100  27206 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   3378563.4   933758.3   3.618 0.000421 ***
## time(ts.P1)     -1583.8      464.1  -3.413 0.000855 ***
## mes.February   -14344.4     7785.4  -1.842 0.067650 .  
## mes.March        2947.5     7784.5   0.379 0.705569    
## mes.April       -3832.4     7783.8  -0.492 0.623291    
## mes.May          2088.7     7783.4   0.268 0.788843    
## mes.June        -2509.7     7783.1  -0.322 0.747621    
## mes.July         3645.0     7783.0   0.468 0.640318    
## mes.August       3959.5     7783.1   0.509 0.611792    
## mes.September   -5034.1     7783.4  -0.647 0.518903    
## mes.October      3042.1     7783.8   0.391 0.696561    
## mes.November    -3062.5     7784.5  -0.393 0.694648    
## mes.December     3977.4     7785.4   0.511 0.610286    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 19440 on 132 degrees of freedom
## Multiple R-squared:  0.1386, Adjusted R-squared:  0.06033 
## F-statistic:  1.77 on 12 and 132 DF,  p-value: 0.05938

mes. <- season(ts.Q2)
reg.Q <- lm(ts.Q2 ~ time(ts.Q2) + mes.)
summary(reg.Q)

## 
## Call:
## lm(formula = ts.Q2 ~ time(ts.Q2) + mes.)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1690.66  -434.73    63.37   513.65  1103.79 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   632777.52   30969.83  20.432  < 2e-16 ***
## time(ts.Q2)     -308.97      15.39 -20.073  < 2e-16 ***
## mes.February    -920.38     258.22  -3.564 0.000508 ***
## mes.March        170.02     258.19   0.659 0.511359    
## mes.April       -213.37     258.17  -0.826 0.410028    
## mes.May          162.77     258.15   0.631 0.529438    
## mes.June        -123.94     258.14  -0.480 0.631931    
## mes.July         198.84     258.14   0.770 0.442511    
## mes.August       111.50     258.14   0.432 0.666489    
## mes.September   -485.78     258.15  -1.882 0.062070 .  
## mes.October     -334.29     258.17  -1.295 0.197632    
## mes.November    -789.18     258.19  -3.057 0.002710 ** 
## mes.December    -281.16     258.22  -1.089 0.278208    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 644.8 on 132 degrees of freedom
## Multiple R-squared:  0.7742, Adjusted R-squared:  0.7537 
## F-statistic: 37.73 on 12 and 132 DF,  p-value: < 2.2e-16

Aplicamos las transformaciones necesarios (ej. diferencias y log)

ts.ddP1 <- diff(diff(ts.P1,12)) # diferencia de produccion
ts.dQ2 <- diff(diff(ts.Q2,12)) # diferencia de produccion

ggtsdisplay(ts.ddP1)

ggtsdisplay(ts.dQ2)

summary(ur.df(ts.ddP1))

## 
## ############################################### 
## # 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 
## -19581.8  -2982.5   -215.1   2861.5  15732.0 
## 
## Coefficients:
##            Estimate Std. Error t value Pr(>|t|)    
## z.lag.1    -1.40162    0.13202 -10.617   <2e-16 ***
## z.diff.lag  0.21329    0.08679   2.458   0.0153 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5477 on 128 degrees of freedom
## Multiple R-squared:  0.5952, Adjusted R-squared:  0.5888 
## F-statistic: 94.09 on 2 and 128 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic is: -10.6165 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau1 -2.58 -1.95 -1.62

summary(ur.df(ts.dQ2))

## 
## ############################################### 
## # 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 
## -1481.95  -266.51   -47.95   278.47  1172.65 
## 
## Coefficients:
##            Estimate Std. Error t value Pr(>|t|)    
## z.lag.1    -1.50675    0.14185 -10.622   <2e-16 ***
## z.diff.lag  0.14661    0.08761   1.673   0.0967 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 447 on 128 degrees of freedom
## Multiple R-squared:  0.6635, Adjusted R-squared:  0.6583 
## F-statistic: 126.2 on 2 and 128 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic is: -10.622 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau1 -2.58 -1.95 -1.62

Se deben guardar ambas variables en un único objeto.

base3 <- cbind.zoo(p = ts.ddP1, q = ts.dQ2) # ocupamos función 'zoo' del paquete del mismo nombre.
View(base3) # la serie está cortada por las diferencias
autoplot(base3) # la serie está cortada por las diferencias

Las variables tienen la misma longitud, por lo tanto no restringimos la muestra.

Ocupamos ‘base3’ para comenzar nuestro análisis VAR.

VARselect(base3, lag.max = 24, type = "const")[["selection"]]

## AIC(n)  HQ(n)  SC(n) FPE(n) 
##     12      1      1     12

VARselect(base3, lag.max = 24, type = "const")

## $selection
## AIC(n)  HQ(n)  SC(n) FPE(n) 
##     12      1      1     12 
## 
## $criteria
##                   1            2            3            4            5
## AIC(n) 2.936836e+01 2.932897e+01 2.936392e+01 2.939296e+01 2.944193e+01
## HQ(n)  2.942877e+01 2.942967e+01 2.950489e+01 2.957421e+01 2.966346e+01
## SC(n)  2.951736e+01 2.957732e+01 2.971160e+01 2.983998e+01 2.998829e+01
## FPE(n) 5.682351e+12 5.463459e+12 5.659076e+12 5.828220e+12 6.124673e+12
##                   6            7            8            9           10
## AIC(n) 2.946248e+01 2.942400e+01 2.947905e+01 2.949353e+01 2.954309e+01
## HQ(n)  2.972429e+01 2.972608e+01 2.982142e+01 2.987617e+01 2.996601e+01
## SC(n)  3.010818e+01 3.016903e+01 3.032343e+01 3.043724e+01 3.058614e+01
## FPE(n) 6.257672e+12 6.029045e+12 6.380918e+12 6.487679e+12 6.835301e+12
##                  11           12           13           14           15
## AIC(n) 2.952432e+01 2.929593e+01 2.930767e+01 2.936597e+01 2.942727e+01
## HQ(n)  2.998751e+01 2.979941e+01 2.985143e+01 2.995000e+01 3.005158e+01
## SC(n)  3.066671e+01 3.053766e+01 3.064874e+01 3.080637e+01 3.096701e+01
## FPE(n) 6.729731e+12 5.376312e+12 5.464665e+12 5.823710e+12 6.230380e+12
##                  16           17           18           19           20
## AIC(n) 2.946969e+01 2.950785e+01 2.953367e+01 2.957518e+01 2.961227e+01
## HQ(n)  3.013427e+01 3.021271e+01 3.027881e+01 3.036060e+01 3.043797e+01
## SC(n)  3.110877e+01 3.124626e+01 3.137142e+01 3.151227e+01 3.164870e+01
## FPE(n) 6.546870e+12 6.857192e+12 7.102038e+12 7.481001e+12 7.855696e+12
##                  21           22           23           24
## AIC(n) 2.957038e+01 2.959824e+01 2.955400e+01 2.935597e+01
## HQ(n)  3.043636e+01 3.050450e+01 3.050054e+01 3.034278e+01
## SC(n)  3.170615e+01 3.183335e+01 3.188845e+01 3.178976e+01
## FPE(n) 7.633527e+12 7.965839e+12 7.747289e+12 6.472436e+12

var1 <- VAR(base3, p=1, type = "both")
summary(var1)

## 
## VAR Estimation Results:
## ========================= 
## Endogenous variables: p, q 
## Deterministic variables: both 
## Sample size: 131 
## Log Likelihood: -2296.131 
## Roots of the characteristic polynomial:
## 0.3393 0.08881
## Call:
## VAR(y = base3, p = 1, type = "both")
## 
## 
## Estimation results for equation p: 
## ================================== 
## p = p.l1 + q.l1 + const + trend 
## 
##         Estimate Std. Error t value Pr(>|t|)
## p.l1    -0.13908    0.09002  -1.545    0.125
## q.l1    -0.95721    1.07298  -0.892    0.374
## const -526.54117  997.46790  -0.528    0.599
## trend    3.26651   12.96374   0.252    0.801
## 
## 
## Residual standard error: 5602 on 127 degrees of freedom
## Multiple R-Squared: 0.03134, Adjusted R-squared: 0.008462 
## F-statistic:  1.37 on 3 and 127 DF,  p-value: 0.255 
## 
## 
## Estimation results for equation q: 
## ================================== 
## q = p.l1 + q.l1 + const + trend 
## 
##        Estimate Std. Error t value Pr(>|t|)   
## p.l1  -0.010515   0.007239  -1.453  0.14883   
## q.l1  -0.289011   0.086284  -3.350  0.00107 **
## const 34.179076  80.211987   0.426  0.67075   
## trend -0.635686   1.042487  -0.610  0.54310   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 450.5 on 127 degrees of freedom
## Multiple R-Squared: 0.1158,  Adjusted R-squared: 0.09494 
## F-statistic: 5.545 on 3 and 127 DF,  p-value: 0.00131 
## 
## 
## 
## Covariance matrix of residuals:
##          p      q
## p 31378152 502442
## q   502442 202912
## 
## Correlation matrix of residuals:
##        p      q
## p 1.0000 0.1991
## q 0.1991 1.0000

causality(var1, cause='p')$Granger

## 
##  Granger causality H0: p do not Granger-cause q
## 
## data:  VAR object var1
## F-Test = 2.1098, df1 = 1, df2 = 254, p-value = 0.1476

causality(var1, cause='q')$Granger

## 
##  Granger causality H0: q do not Granger-cause p
## 
## data:  VAR object var1
## F-Test = 0.79585, df1 = 1, df2 = 254, p-value = 0.3732

plot(irf(var1, impulse='p', response = 'q'))

plot(irf(var1, impulse='q', response = 'p'))

var.serial <- serial.test(var1, lags.pt = 12)

var.serial

## 
##  Portmanteau Test (asymptotic)
## 
## data:  Residuals of VAR object var1
## Chi-squared = 121.36, df = 44, p-value = 3.651e-09

plot(var.serial)

Otra propuesta

VARselect(base3, lag.max = 24, type = "const")[["selection"]]

## AIC(n)  HQ(n)  SC(n) FPE(n) 
##     12      1      1     12

var12 <- VAR(base3, p=12, type = "both")
summary(var12)

## 
## VAR Estimation Results:
## ========================= 
## Endogenous variables: p, q 
## Deterministic variables: both 
## Sample size: 120 
## Log Likelihood: -2042.279 
## Roots of the characteristic polynomial:
## 0.9827 0.9827 0.9689 0.9689 0.9579 0.9579 0.9478 0.9478 0.9381 0.9381 0.9377 0.9377 0.9321 0.9321 0.9153 0.9153 0.9001 0.9001 0.8875 0.8875 0.8577 0.8577 0.7322 0.7322
## Call:
## VAR(y = base3, p = 12, type = "both")
## 
## 
## Estimation results for equation p: 
## ================================== 
## p = p.l1 + q.l1 + p.l2 + q.l2 + p.l3 + q.l3 + p.l4 + q.l4 + p.l5 + q.l5 + p.l6 + q.l6 + p.l7 + q.l7 + p.l8 + q.l8 + p.l9 + q.l9 + p.l10 + q.l10 + p.l11 + q.l11 + p.l12 + q.l12 + const + trend 
## 
##         Estimate Std. Error t value Pr(>|t|)   
## p.l1  -1.291e-01  1.011e-01  -1.276   0.2051   
## q.l1  -1.608e+00  1.298e+00  -1.239   0.2183   
## p.l2  -1.482e-01  1.026e-01  -1.445   0.1519   
## q.l2  -1.145e+00  1.353e+00  -0.846   0.3997   
## p.l3  -9.188e-03  1.041e-01  -0.088   0.9299   
## q.l3  -2.034e-01  1.354e+00  -0.150   0.8809   
## p.l4  -2.454e-02  1.063e-01  -0.231   0.8180   
## q.l4  -1.212e+00  1.358e+00  -0.892   0.3746   
## p.l5  -1.022e-01  1.139e-01  -0.897   0.3720   
## q.l5  -1.738e+00  1.376e+00  -1.263   0.2097   
## p.l6  -3.300e-02  1.120e-01  -0.295   0.7688   
## q.l6  -2.730e+00  1.409e+00  -1.937   0.0557 . 
## p.l7   1.435e-01  1.110e-01   1.293   0.1993   
## q.l7  -2.645e+00  1.383e+00  -1.912   0.0589 . 
## p.l8  -3.282e-02  1.101e-01  -0.298   0.7663   
## q.l8   1.178e-01  1.416e+00   0.083   0.9339   
## p.l9   4.893e-02  1.101e-01   0.445   0.6577   
## q.l9   1.142e+00  1.393e+00   0.820   0.4142   
## p.l10  2.147e-02  1.095e-01   0.196   0.8449   
## q.l10  2.113e-01  1.369e+00   0.154   0.8777   
## p.l11  8.771e-02  1.073e-01   0.817   0.4159   
## q.l11  8.827e-01  1.368e+00   0.645   0.5202   
## p.l12 -2.996e-01  1.084e-01  -2.762   0.0069 **
## q.l12  1.109e+00  1.282e+00   0.865   0.3891   
## const -4.426e+02  1.125e+03  -0.393   0.6949   
## trend -8.224e-01  1.402e+01  -0.059   0.9534   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 5076 on 94 degrees of freedom
## Multiple R-Squared: 0.3157,  Adjusted R-squared: 0.1337 
## F-statistic: 1.735 on 25 and 94 DF,  p-value: 0.0306 
## 
## 
## Estimation results for equation q: 
## ================================== 
## q = p.l1 + q.l1 + p.l2 + q.l2 + p.l3 + q.l3 + p.l4 + q.l4 + p.l5 + q.l5 + p.l6 + q.l6 + p.l7 + q.l7 + p.l8 + q.l8 + p.l9 + q.l9 + p.l10 + q.l10 + p.l11 + q.l11 + p.l12 + q.l12 + const + trend 
## 
##         Estimate Std. Error t value Pr(>|t|)    
## p.l1  -0.0118261  0.0075501  -1.566  0.12063    
## q.l1  -0.2999644  0.0968522  -3.097  0.00258 ** 
## p.l2   0.0070491  0.0076597   0.920  0.35978    
## q.l2  -0.1462749  0.1010088  -1.448  0.15091    
## p.l3  -0.0038680  0.0077739  -0.498  0.61996    
## q.l3  -0.1165660  0.1010569  -1.153  0.25164    
## p.l4  -0.0119375  0.0079370  -1.504  0.13593    
## q.l4  -0.0626822  0.1013965  -0.618  0.53795    
## p.l5   0.0002905  0.0085001   0.034  0.97281    
## q.l5  -0.0333277  0.1027013  -0.325  0.74627    
## p.l6  -0.0088207  0.0083567  -1.056  0.29389    
## q.l6  -0.0157291  0.1051925  -0.150  0.88146    
## p.l7   0.0080932  0.0082854   0.977  0.33118    
## q.l7   0.0378268  0.1032574   0.366  0.71494    
## p.l8  -0.0041707  0.0082166  -0.508  0.61293    
## q.l8   0.0232301  0.1056707   0.220  0.82648    
## p.l9  -0.0006136  0.0082149  -0.075  0.94062    
## q.l9  -0.0011511  0.1039794  -0.011  0.99119    
## p.l10 -0.0041929  0.0081706  -0.513  0.60903    
## q.l10  0.1439504  0.1021882   1.409  0.16223    
## p.l11  0.0098111  0.0080120   1.225  0.22381    
## q.l11  0.0767256  0.1020928   0.752  0.45421    
## p.l12  0.0118739  0.0080948   1.467  0.14575    
## q.l12 -0.3976084  0.0957015  -4.155 7.17e-05 ***
## const 91.1283805 83.9899431   1.085  0.28070    
## trend -1.3586284  1.0468408  -1.298  0.19752    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 378.9 on 94 degrees of freedom
## Multiple R-Squared: 0.4814,  Adjusted R-squared: 0.3434 
## F-statistic:  3.49 on 25 and 94 DF,  p-value: 6.005e-06 
## 
## 
## 
## Covariance matrix of residuals:
##          p      q
## p 25763585 558880
## q   558880 143544
## 
## Correlation matrix of residuals:
##        p      q
## p 1.0000 0.2906
## q 0.2906 1.0000

causality(var12, cause='p')$Granger

## 
##  Granger causality H0: p do not Granger-cause q
## 
## data:  VAR object var12
## F-Test = 0.90703, df1 = 12, df2 = 188, p-value = 0.541

causality(var12, cause='q')$Granger

## 
##  Granger causality H0: q do not Granger-cause p
## 
## data:  VAR object var12
## F-Test = 1.2406, df1 = 12, df2 = 188, p-value = 0.2581

plot(irf(var12, impulse='p', response = 'q'))

plot(irf(var12, impulse='q', response = 'p'))

var.serial <- serial.test(var12, lags.pt = 12)
var.serial

## 
##  Portmanteau Test (asymptotic)
## 
## data:  Residuals of VAR object var12
## Chi-squared = 18.823, df = 0, p-value < 2.2e-16

plot(var.serial)

14. Evaluar la potencial cointegracion usando la metodologia de Engle-Granger y Johan-Sen

rm(list = ls())
base2<-read_xls("C:/Users/rodar/OneDrive/Escritorio/SERIES/PNG Y LIQ.xls")

ts.P<-ts(base2$P, start = c(2006,1), frequency = 12)
ts.Q<-ts(base2$Q, start = c(2006,1), frequency = 12)
par(mfrow=c(1,2))
plot(ts.P, main = "PRODUCCION NACIONAL GAS")
plot(ts.Q, main = "PRODUCCION NACIONAL METANO")

ur.df(ts.Q)

## 
## ############################################################### 
## # Augmented Dickey-Fuller Test Unit Root / Cointegration Test # 
## ############################################################### 
## 
## The value of the test statistic is: -1.3179

ur.df(ts.P)

## 
## ############################################################### 
## # Augmented Dickey-Fuller Test Unit Root / Cointegration Test # 
## ############################################################### 
## 
## The value of the test statistic is: -0.1038

Evidencia de raiz unitaria en ambas variables.

lr.reg<-dynlm(ts.Q~L(ts.Q)+ ts.P+L(ts.Q))
summary(lr.reg)

## 
## Time series regression with "ts" data:
## Start = 2006(2), End = 2018(1)
## 
## Call:
## dynlm(formula = ts.Q ~ L(ts.Q) + ts.P + L(ts.Q))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1260.04  -382.88    -6.49   297.08  1767.79 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -3.051e+02  5.661e+02  -0.539    0.591    
## L(ts.Q)      8.370e-01  3.965e-02  21.108  < 2e-16 ***
## ts.P         1.070e-02  2.549e-03   4.197 4.77e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 586.9 on 141 degrees of freedom
## Multiple R-squared:  0.7942, Adjusted R-squared:  0.7913 
## F-statistic: 272.1 on 2 and 141 DF,  p-value: < 2.2e-16

dwtest(lr.reg)

## 
##  Durbin-Watson test
## 
## data:  lr.reg
## DW = 2.607, p-value = 0.9998
## alternative hypothesis: true autocorrelation is greater than 0

Adicionalmente no hay correlacion serial.

error<-residuals(lr.reg)
plot(error)
abline(h=0)

ur.df(error)

## 
## ############################################################### 
## # Augmented Dickey-Fuller Test Unit Root / Cointegration Test # 
## ############################################################### 
## 
## The value of the test statistic is: -8.156

El valor critico al 5% es de -8.156, por lo que no hay raiz unitaria

datar<- ts.intersect(ts.Q, ts.P, error)
dyn<- dynlm(d(ts.Q)~L(error)+L(d(ts.P))+L(d(ts.Q)), data = datar)
summary(dyn)

## 
## Time series regression with "ts" data:
## Start = 2006(4), End = 2018(1)
## 
## Call:
## dynlm(formula = d(ts.Q) ~ L(error) + L(d(ts.P)) + L(d(ts.Q)), 
##     data = datar)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1348.85  -281.08    55.86   340.11  1260.65 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)  
## (Intercept) -43.663604  43.821386  -0.996   0.3208  
## L(error)     -0.102393   0.190340  -0.538   0.5915  
## L(d(ts.P))   -0.020216   0.007752  -2.608   0.0101 *
## L(d(ts.Q))   -0.190941   0.201984  -0.945   0.3461  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 515.5 on 138 degrees of freedom
## Multiple R-squared:  0.3306, Adjusted R-squared:  0.316 
## F-statistic: 22.72 on 3 and 138 DF,  p-value: 5.159e-12

El procedimiento de Engle y Granger puede resumirse en la prueba de Phillips-oularis.

library(tseries)
po.test(datar[,1:2])

## Warning in po.test(datar[, 1:2]): p-value greater than printed p-value

## 
##  Phillips-Ouliaris Cointegration Test
## 
## data:  datar[, 1:2]
## Phillips-Ouliaris demeaned = -8.1421, Truncation lag parameter =
## 1, p-value = 0.15

ANALISIS ARIMA DE LA PRODUCCION DE GAS