Indicaciones: El taller debe ser desarrollado en R. Una vez terminado,debe ser disponibilizado en la plataforma.

Punto 1.

Explore los datos del paquete . Obtenga los soportes visuales y los aspectos descriptivos de la serie a partir de los comandos dados.

#install.packages("astsa")
library(astsa)
data(birth)
birth
##      Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 1948 295 286 300 278 272 268 308 321 313 308 291 296
## 1949 294 273 300 271 282 285 318 323 313 311 291 293
## 1950 297 273 294 259 276 294 316 325 315 312 292 301
## 1951 304 282 313 296 313 307 328 334 329 329 304 312
## 1952 312 300 317 292 300 311 345 350 344 336 315 323
## 1953 322 296 315 287 307 321 354 356 348 334 320 340
## 1954 332 302 324 305 318 329 359 363 359 352 335 342
## 1955 329 306 332 309 326 325 354 367 362 354 337 345
## 1956 339 325 345 309 315 334 370 383 375 370 344 355
## 1957 346 317 348 331 345 348 380 381 377 376 348 356
## 1958 344 320 347 326 343 338 361 368 378 374 347 358
## 1959 349 323 358 331 338 343 374 380 377 368 346 358
## 1960 338 329 347 327 335 336 370 399 385 368 351 362
## 1961 358 333 356 335 348 346 374 386 384 372 343 346
## 1962 346 318 359 328 333 329 366 373 367 363 337 346
## 1963 355 314 343 322 336 327 362 366 361 358 327 330
## 1964 336 326 337 316 331 331 359 350 356 347 328 336
## 1965 315 292 322 291 302 310 330 335 333 318 305 313
## 1966 301 281 302 291 297 291 311 319 317 317 296 307
## 1967 295 265 300 271 291 290 310 318 310 304 285 288
## 1968 277 260 282 274 288 287 308 312 306 304 282 305
## 1969 284 273 286 284 294 288 315 322 317 309 295 306
## 1970 300 275 301 292 298 306 326 332 329 328 308 324
## 1971 299 284 306 290 292 285 295 306 317 305 294 287
## 1972 278 261 275 256 270 264 265 284 284 275 269 275
## 1973 259 244 267 255 260 253 267 277 277 264 255 260
## 1974 261 238 257 246 254 255 273 276 286 283 261 276
## 1975 264 243 259 250 262 253 280 288 270 273 241 266
## 1976 257 242 266 241 252 250 281 278 286 278 260 272
## 1977 274 256 276 259 273 272 297 296 290 282 262 275
## 1978 262 251 285 260 272 265 296 312 289 282 274 281
## 1979 277
library(TSstudio)
## Warning: package 'TSstudio' was built under R version 4.1.2
library(highcharter)
## Warning: package 'highcharter' was built under R version 4.1.2
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(foreign)
library(ggfortify)
## Warning: package 'ggfortify' was built under R version 4.1.2
## Loading required package: ggplot2
library(forecast)
## Warning: package 'forecast' was built under R version 4.1.2
## Registered S3 methods overwritten by 'forecast':
##   method                 from     
##   autoplot.Arima         ggfortify
##   autoplot.acf           ggfortify
##   autoplot.ar            ggfortify
##   autoplot.bats          ggfortify
##   autoplot.decomposed.ts ggfortify
##   autoplot.ets           ggfortify
##   autoplot.forecast      ggfortify
##   autoplot.stl           ggfortify
##   autoplot.ts            ggfortify
##   fitted.ar              ggfortify
##   fortify.ts             ggfortify
##   residuals.ar           ggfortify
## 
## Attaching package: 'forecast'
## The following object is masked from 'package:astsa':
## 
##     gas
library(gridExtra)
## Warning: package 'gridExtra' was built under R version 4.1.2
library(seasonal)
## Warning: package 'seasonal' was built under R version 4.1.2
## 
## Attaching package: 'seasonal'
## The following object is masked from 'package:astsa':
## 
##     unemp
library(lattice)
library(zoo)
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(urca)
## Warning: package 'urca' was built under R version 4.1.2
library(dynlm)
## Warning: package 'dynlm' was built under R version 4.1.2
library(ggplot2)
ts_decompose(birth)
autoplot(birth,ts.colour = "red")+
  ggtitle("Nacimientos mensuales 1948 - 1979")+
  xlab("Tiempo")+
  ylab("Número de nacimientos")
## Warning: Ignoring unknown parameters: ts.colour

plot(birth)

ts_plot(birth)
autoplot(birth)

hchart(birth)
hchart(stl(birth, s.window='periodic'))
## Warning: Deprecated function. Use the `create_axis` function.
ts_seasonal(birth, type = "all")
monthplot(birth)

Descripción de los datos Punto 1

La base esta compuesta por datos mensuales desde Enero 1948 hasta enero 1979

Tendencia: para la serie de nacimientos graficamente no se ve una tendencia pues los primeros años de la serie se ve un crecimiento hasta 1960 donde inicia un decrecimiento de los nacimiento, crece nuevamente en 1978

Estacionalidad: la serie es mensual vemos un comportamiento estacional en los meses donde julio, agosto y septiembre son los de mayor número de nacimientos, esto tambien se puede observar en el boxplot donde la distribución de nacimientos mensuales de los años es similar

Punto 2.

Cargue y Explore los datos IPC (archivo de excel adjunto).

Cargue en R la serie con, Fecha de inicio: Enero 2000, Fecha de Finalización: Diciembre de 2019, frecuencia mensual

Obtenga los soportes visuales y los aspectos descriptivos de la serie, obtenga los ciclos de la serie realice un boxplot, y haga un breve análisis.

hchart(stl(IPC, s.window='periodic'))
## Warning: Deprecated function. Use the `create_axis` function.
ts_seasonal(IPC, type = "all")
monthplot(IPC)

start(IPC)
## [1] 2000    1
end(IPC)
## [1] 2019   12
frequency(IPC)
## [1] 12

Descripción y Análisis base IPC

Tendencia: para la serie de IPC graficamente no se ve una tendencia pues vemos la serie constante respecto a la media desde el año 2000 a 2019

Estacionalidad: la serie es mensual si vemos un comportamiento estacional en los meses de enero a febrero con un incremento del IPC explicado por los aumentos de salarios anuales de los trabajadores que se traduce en incrementos de precios de los productos y servios, despues del segundo semestre se observa un decrecimiento del Indide de precios del consumidor hasta el mes de diciembre donde nuevamente inicia el crecimiento

Punto 3

Simule y obtenga los soportes visuales de un Ruido blanco de media 0 realice el acf y pacf. Describa brevemente lo que puede observar.

ruido_blanco=rnorm(1000,0,1)

plot.ts(ruido_blanco, main="Simulación de Ruido Blanco", xlab="Tiempo", ylab="Valores",col="13")

## Media 0
mean(ruido_blanco)
## [1] 0.04379612
## ACF
acf(ruido_blanco)

## PACF
pacf(ruido_blanco)

Descripción y Análisis ruido blanco

para el caso de las series de ruido blanco como son impredecibles no son objeto para poder realizar un modelo, ya detallando el ACF podemos ver que todas las lineas estan dentro de los limites, por lo que no aplicaria un modelo de medias moviles (MA) y en el PAFC se ve que todas las lineas estan dentro de los limites por lo que tampoco se utilizaría un modelo AR

Punto 4.

En sus propias palabras, explique los modelos autorregresivos y genere 100 datos de un modelo cuyos parámetros sean 0.2 y 0.6.

Modelos atoregresivos (AR):modelo estadistico que a partir de datos pasados que conocemos de una serie tiempo, y que la variable de interes depende linealmente de sus valores anteriores podemos precidir comportamientos futuros

layout(matrix(c(1,1,2,3),2,2,byrow = TRUE))
set.seed(123)
ar1 <- arima.sim(n = 100, list(order= c(2,0,0),ar = c(0.2,0.6)))
plot(ar1, main=(expression(AR(2)~~~~theta==c(.2,.6))))
acf(ar1, main="Autocorrelación simple de orden 2")
pacf(ar1, main="Autocorrelación parcial")

\[y_t = 0.2 * y_{t-1}+0.6 * y_{t-2}+ w_t\]

Punto 5

En sus propias palabras, explique los modelos medias móviles y genere 100 datos de un modelo cuyos parámetros sean -0.1 y 0.3.

Modelos de media movil (MA):modelo estadístico que explica el valor de una variable en un periodo t en función de un termino independiente y una sucesión de terminos de error, estos errores corresponden a periodos (t-n) ponderados

layout(matrix(c(1,1,2,3),2,2,byrow = TRUE))
set.seed(123)
ma1 <- arima.sim(n = 100, list(order= c(0,0,2),ma = c(-0.1,0.3)))
plot(ma1, main=(expression(MA(2)~~~~theta==c(-.1,.3))))
acf(ma1, main="Autocorrelación simple de orden 2")
pacf(ma1, main="Autocorrelación parcial")

- Escriba la ecuación del modelo

\[x_t = \mu + w_t +(-0.1)*w_{t-1}+0.3*w_{t-2}\]

Punto 6.

Analizar las siguientes series de tiempo:

y responder a las preguntas:

Paquete Fitar

library(FitAR)
## Warning: package 'FitAR' was built under R version 4.1.2
## Loading required package: leaps
## Warning: package 'leaps' was built under R version 4.1.2
## Loading required package: ltsa
## Loading required package: bestglm
## Warning: package 'bestglm' was built under R version 4.1.2
## 
## Attaching package: 'FitAR'
## The following object is masked from 'package:forecast':
## 
##     BoxCox
data(SeriesA)
SeriesA
## Time Series:
## Start = 1 
## End = 197 
## Frequency = 1 
##   [1] 17.0 16.6 16.3 16.1 17.1 16.9 16.8 17.4 17.1 17.0 16.7 17.4 17.2 17.4 17.4
##  [16] 17.0 17.3 17.2 17.4 16.8 17.1 17.4 17.4 17.5 17.4 17.6 17.4 17.3 17.0 17.8
##  [31] 17.5 18.1 17.5 17.4 17.4 17.1 17.6 17.7 17.4 17.8 17.6 17.5 16.5 17.8 17.3
##  [46] 17.3 17.1 17.4 16.9 17.3 17.6 16.9 16.7 16.8 16.8 17.2 16.8 17.6 17.2 16.6
##  [61] 17.1 16.9 16.6 18.0 17.2 17.3 17.0 16.9 17.3 16.8 17.3 17.4 17.7 16.8 16.9
##  [76] 17.0 16.9 17.0 16.6 16.7 16.8 16.7 16.4 16.5 16.4 16.6 16.5 16.7 16.4 16.4
##  [91] 16.2 16.4 16.3 16.4 17.0 16.9 17.1 17.1 16.7 16.9 16.5 17.2 16.4 17.0 17.0
## [106] 16.7 16.2 16.6 16.9 16.5 16.6 16.6 17.0 17.1 17.1 16.7 16.8 16.3 16.6 16.8
## [121] 16.9 17.1 16.8 17.0 17.2 17.3 17.2 17.3 17.2 17.2 17.5 16.9 16.9 16.9 17.0
## [136] 16.5 16.7 16.8 16.7 16.7 16.6 16.5 17.0 16.7 16.7 16.9 17.4 17.1 17.0 16.8
## [151] 17.2 17.2 17.4 17.2 16.9 16.8 17.0 17.4 17.2 17.2 17.1 17.1 17.1 17.4 17.2
## [166] 16.9 16.9 17.0 16.7 16.9 17.3 17.8 17.8 17.6 17.5 17.0 16.9 17.1 17.2 17.4
## [181] 17.5 17.9 17.0 17.0 17.0 17.2 17.3 17.4 17.4 17.0 18.0 18.2 17.6 17.8 17.7
## [196] 17.2 17.4
## attr(,"title")
## [1] Chemical process concentrations
library(car)
## Warning: package 'car' was built under R version 4.1.2
## Loading required package: carData
## Warning: package 'carData' was built under R version 4.1.2
## 
## Attaching package: 'car'
## The following object is masked from 'package:FitAR':
## 
##     Boot
library(urca)
library(forecast)
library(tseries)
## Warning: package 'tseries' was built under R version 4.1.2
adf.test(SeriesA)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  SeriesA
## Dickey-Fuller = -2.6562, Lag order = 5, p-value = 0.3014
## alternative hypothesis: stationary

Aplicando la prueba de Dickey-Fuller aumentado obtenemos un p-valor de 0.30 que es mayor al 0.05 por lo tanto rechazamos la Ho por lo que la serie no es estacionaria por lo que debemos aplicarle una diferencia

adf.test(diff(SeriesA))
## Warning in adf.test(diff(SeriesA)): p-value smaller than printed p-value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  diff(SeriesA)
## Dickey-Fuller = -9.9271, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationary
par(mfrow=c(2,3))  # garaficar en 2 filas y tres columnas
plot(SeriesA,type="o")
acf(SeriesA)
pacf(SeriesA)
plot(diff(SeriesA),type="o")
abline(h=2*sqrt(var(diff(SeriesA))),col="red",lty=2)
abline(h=-2*sqrt(var(diff(SeriesA))),col="red",lty=2)
acf(diff(SeriesA))
pacf(diff(SeriesA))

Con el resultado vamos a probar tres modelos ARIMA (3,1,1), ARIMA (1,1,1) y un ARIMA (2,1,1)

Modelo1: ARIMA (3,1,1)

modelo1<-stats::arima(SeriesA,
                      order=c(3,1,1), fixed=c(NA,NA,NA,NA))
modelo1
## 
## Call:
## stats::arima(x = SeriesA, order = c(3, 1, 1), fixed = c(NA, NA, NA, NA))
## 
## Coefficients:
##          ar1     ar2      ar3      ma1
##       0.2201  0.0834  -0.0582  -0.8331
## s.e.  0.1118  0.0926   0.0864   0.0885
## 
## sigma^2 estimated as 0.0977:  log likelihood = -50.57,  aic = 111.15
tt <- modelo1$coef[which(modelo1$coef!=0)]/sqrt(diag(modelo1$var.coef))
tt
##        ar1        ar2        ar3        ma1 
##  1.9687450  0.9015403 -0.6732022 -9.4166923
1 - pt(abs(tt),(modelo1$nobs - length(modelo1$coef[which(modelo1$coef!=0)])))
##       ar1       ar2       ar3       ma1 
## 0.0252107 0.1842153 0.2508141 0.0000000
BIC (modelo1)
## [1] 127.5375

Modelo2: ARIMA (1,1,1)

modelo2<-stats::arima(SeriesA,
                      order=c(1,1,1), fixed=c(NA,NA))
modelo2
## 
## Call:
## stats::arima(x = SeriesA, order = c(1, 1, 1), fixed = c(NA, NA))
## 
## Coefficients:
##          ar1      ma1
##       0.2155  -0.8193
## s.e.  0.1011   0.0631
## 
## sigma^2 estimated as 0.09851:  log likelihood = -51.37,  aic = 108.74
tt2 <- modelo2$coef[which(modelo2$coef!=0)]/sqrt(diag(modelo2$var.coef))
tt2
##        ar1        ma1 
##   2.130938 -12.993605
1 - pt(abs(tt2),(modelo2$nobs - length(modelo2$coef[which(modelo2$coef!=0)])))
##        ar1        ma1 
## 0.01717663 0.00000000
BIC (modelo2)
## [1] 118.5766

Modelo3: ARIMA (2,1,1)

modelo3<-stats::arima(SeriesA,
                      order=c(2,1,1), fixed=c(NA,NA,NA))
modelo3
## 
## Call:
## stats::arima(x = SeriesA, order = c(2, 1, 1), fixed = c(NA, NA, NA))
## 
## Coefficients:
##          ar1     ar2      ma1
##       0.2469  0.0946  -0.8639
## s.e.  0.0978  0.0881   0.0636
## 
## sigma^2 estimated as 0.09792:  log likelihood = -50.8,  aic = 109.6
tt3 <- modelo3$coef[which(modelo3$coef!=0)]/sqrt(diag(modelo3$var.coef))
tt3
##        ar1        ar2        ma1 
##   2.523817   1.073214 -13.576199
1 - pt(abs(tt3),(modelo3$nobs - length(modelo3$coef[which(modelo3$coef!=0)])))
##         ar1         ar2         ma1 
## 0.006207457 0.142257504 0.000000000
BIC (modelo3)
## [1] 122.7094

Los BIC de los modelos son:

El mejor modelo es el ARIMA (1,1,1) por tener el menor BIC

Ecuación

\[x_t = \mu + \triangle^{1}w_t +(-12,9)*\triangle^{1}w_{t-1}+2.13 * \triangle^{1}y_{t-1}\]

plot(forecast(modelo2,h=4, fan=T))
lines(fitted(modelo2), col="red")

hchart(forecast(modelo2,h=4, fan=T))

Descripción del fenomeno: la linea roja que muestra que el comportamiento del modelo 2 se asemeja al comportamiento de la serie original y los 4 datos que proyectamos tienen un crecimiento leve pero esta entre los limites de confianza

Datos de blackboard

library(readxl)
data <- read.table("C:/Users/juanc/Downloads/data (1).txt", header = TRUE)
View(data)
dt= ts(data,start = c(1604,1), frequency = 12)
adf.test(dt)
## Warning in adf.test(dt): p-value smaller than printed p-value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  dt
## Dickey-Fuller = -16.95, Lag order = 17, p-value = 0.01
## alternative hypothesis: stationary

Aplicando la prueba de Dickey-Fuller aumentado obtenemos un p-valor de 0.01 que es menor al 0.05 por lo tanto no rechazamos la Ho por lo que la serie es estacionaria

layout(matrix(c(1,1,2,3),2,2,byrow = TRUE))
set.seed(123)

plot(dt, main="Base Data.txt")
acf(dt, main="Autocorrelación simple")
pacf(dt, main="Autocorrelación parcial")

Con el resultado vamos a probar tres modelos ARIMA (8,0,1), ARIMA (3,0,1) y un ARIMA (4,0,1)

Modelo1: ARIMA (8,0,1)

modelo1<-stats::arima(dt,
                    order=c(8,0,1), fixed=c(NA,NA,NA,NA,NA,NA,NA,NA,NA,NA))
modelo1
## 
## Call:
## stats::arima(x = dt, order = c(8, 0, 1), fixed = c(NA, NA, NA, NA, NA, NA, NA, 
##     NA, NA, NA))
## 
## Coefficients:
##          ar1      ar2      ar3      ar4     ar5      ar6     ar7      ar8
##       0.1080  -0.0082  -0.0114  -0.0023  0.0394  -0.0239  0.0256  -0.0032
## s.e.  0.0402   0.0368   0.0307   0.0253  0.0220   0.0207  0.0192   0.0175
##          ma1  intercept
##       0.7986    -0.0067
## s.e.  0.0377     0.0298
## 
## sigma^2 estimated as 1.052:  log likelihood = -7221.39,  aic = 14464.78
tt <- modelo1$coef[which(modelo1$coef!=0)]/sqrt(diag(modelo1$var.coef))
tt
##         ar1         ar2         ar3         ar4         ar5         ar6 
##  2.68525244 -0.22264473 -0.37104493 -0.09068503  1.79079064 -1.15681505 
##         ar7         ar8         ma1   intercept 
##  1.33685627 -0.18437672 21.18251174 -0.22457656
1 - pt(abs(tt),(modelo1$nobs - length(modelo1$coef[which(modelo1$coef!=0)])))
##         ar1         ar2         ar3         ar4         ar5         ar6 
## 0.003635738 0.411910567 0.355309908 0.463873254 0.036693746 0.123701653 
##         ar7         ar8         ma1   intercept 
## 0.090665250 0.426862711 0.000000000 0.411158948
BIC (modelo1)
## [1] 14536.47

Modelo2: ARIMA (3,0,1)

modelo2<-stats::arima(dt,
                      order=c(3,0,1), fixed=c(NA,NA,NA,NA,NA))
modelo2
## 
## Call:
## stats::arima(x = dt, order = c(3, 0, 1), fixed = c(NA, NA, NA, NA, NA))
## 
## Coefficients:
##          ar1     ar2      ar3     ma1  intercept
##       0.0726  0.0232  -0.0343  0.8334    -0.0068
## s.e.  0.0197  0.0190   0.0171  0.0135     0.0283
## 
## sigma^2 estimated as 1.053:  log likelihood = -7224.82,  aic = 14461.65
tt2 <- modelo2$coef[which(modelo2$coef!=0)]/sqrt(diag(modelo2$var.coef))
tt2
##        ar1        ar2        ar3        ma1  intercept 
##  3.6797463  1.2192754 -2.0018400 61.6233146 -0.2387355
1 - pt(abs(tt2),(modelo2$nobs - length(modelo2$coef[which(modelo2$coef!=0)])))
##          ar1          ar2          ar3          ma1    intercept 
## 0.0001179631 0.1113986184 0.0226779619 0.0000000000 0.4056602218
BIC (modelo2)
## [1] 14500.75

Modelo3: ARIMA (4,0,1)

modelo3<-stats::arima(dt,
                      order=c(4,0,1), fixed=c(NA,NA,NA,NA,NA,NA))
modelo3
## 
## Call:
## stats::arima(x = dt, order = c(4, 0, 1), fixed = c(NA, NA, NA, NA, NA, NA))
## 
## Coefficients:
##          ar1     ar2      ar3     ar4     ma1  intercept
##       0.0612  0.0338  -0.0448  0.0253  0.8455    -0.0067
## s.e.  0.0205  0.0198   0.0182  0.0165  0.0147     0.0290
## 
## sigma^2 estimated as 1.053:  log likelihood = -7223.65,  aic = 14461.3
tt3 <- modelo3$coef[which(modelo3$coef!=0)]/sqrt(diag(modelo3$var.coef))
tt3
##        ar1        ar2        ar3        ar4        ma1  intercept 
##  2.9785789  1.7110985 -2.4612776  1.5320993 57.3515211 -0.2312567
1 - pt(abs(tt3),(modelo3$nobs - length(modelo3$coef[which(modelo3$coef!=0)])))
##         ar1         ar2         ar3         ar4         ma1   intercept 
## 0.001454906 0.043562517 0.006938948 0.062780639 0.000000000 0.408562423
BIC (modelo3)
## [1] 14506.92

Los BIC de los modelos son:

-ARIMA (8,0,1):14536.47 -ARIMA (3,0,1):14500.75 -ARIMA (4,0,1):14506.92

El mejor modelo es el ARIMA (3,0,1) por tener el menor BIC

Ecuación

\[x_t = \mu + w_t +61,6*w_{t-1}+3,7 * y_{t-1}+1,2 * y_{t-2}+(-2) * y_{t-3}\]

plot(forecast(modelo2,h=4, fan=T))
lines(fitted(modelo2), col="red")

hchart(forecast(modelo2,h=4, fan=T))

Descripción del fenomeno: la linea roja que muestra que el comportamiento del modelo 2 se asemeja al comportamiento de la serie original (estacionaria) y los 4 datos que proyectamos tienen un comportamiento similar estando entre los limites de confianza