¿Qué son las series de tiempo según autores?

Las series de tiempo son colecciones de observaciones sobre un determinado fenómeno efectuadas en sucesivos momentos del tiempo, usualmente equiespaciados . Corresponde a una realización de un proceso generador de datos.

¿Cuáles son los 4 componentes de una serie de tiempo?

Estas cuatro componentes son: Tendencia secular, variación estacional, variación cíclica y variación irregular . Supondremos, además, que existe una relación multiplicativa entre estas cuatro componentes; es decir, cualquier valor de una serie es el producto de factores que se pueden atribuir a las cuatro componentes.

¿Cuál es el modelo matemático más utilizado de una serie de tiempo?

El modelo matemático de serie temporal más utilizado es el modelo de media móvil integrada autorregresiva (ARIMA) . Este modelo se utiliza ampliamente en diversos campos, como la economía, las finanzas y la previsión, para analizar y predecir valores futuros basándose en patrones de datos pasados.

[PDF] Notas de Clase Series de Tiempo con R

Download Notas de Clase Series de Tiempo con R...

Notas de Clase Series de Tiempo con R

NOTAS DE CLASE Series de Tiempo con R

´ NORMAN GIRALDO GOMEZ Profesor Asociado Escuela de Estad´ıstica Universidad Nacional de Colombia Medell´ın

Universidad Nacional de Colombia Medellín

c Copyright 2006 Universidad Nacional de Colombia. Profesor Norman Diego Giraldo Gómez. Publicado por Editorial ISBN 000-000-000-0 No está permitido reproducir esta publicación o transmitirla por cualquier forma o medio, electrónico, mecánico, fotocopiado, escaneo ó de otro tipo excepto para citas cortas, sin el permiso de la Editorial.

Centro de Documentación Rafael Botero, UN-Medellín: Series de Tiempo con R/ Norman Diego Giraldo G´omez. p. cm.—(Colecci´on Notas de Clase) “Universidad Nacional de Colombia." Incluye referencias bibliogr´aficas e ´ındice. ISBN 0-000-00000-0 (pbk.) 1. Probabilidades—Teor´ıa. 2. Matem´aticas Ciencias—Investigaci´on—Teor´ıa. I. Giraldo, Norman D. II. Series.

519.2 G887c Diagramación en LaTeX. Impresión Editorial

´Indice general

1. Introduci´on al Lenguaje R

1.1. Introducci´on al Curso . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1.1.1. Notas sobre Pron´osticos . . . . . . . . . . . . . . . . . . . . . . .

1.2. Introducci´on al R . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1.3. Uso de R de Manera Interactiva versus Ejecuci´on desde un Programa . . . .

1.3.1. Interactiva . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1.3.2. Ejecuci´on desde un Programa . . . . . . . . . . . . . . . . . . . .

1.3.3. Lectura de Datos . . . . . . . . . . . . . . . . . . . . . . . . . . .

1.4. Repaso del Modelo de Regresi´on Lineal . . . . . . . . . . . . . . . . . . .

1.4.1. Estimaci´on de los Par´ametros . . . . . . . . . . . . . . . . . . . .

1.4.2. Pruebas de Ajuste del Modelo . . . . . . . . . . . . . . . . . . . .

1.4.3. Estad´ısticos de Ajuste y de Selecci´on de Modelos . . . . . . . . . .

1.4.4. M´ınimos Cuadrados Nolineales . . . . . . . . . . . . . . . . . . .

1.5. Ejemplo de Regresi´on con R . . . . . . . . . . . . . . . . . . . . . . . . .

1.6. Manejo de Series de Tiempo con R . . . . . . . . . . . . . . . . . . . . . .

vi 2. Modelado y Pron´ostico de la Tendencia

2.1. El Modelo Aditivo de Componentes de Series de Tiempo . . . . . . . . . .

2.2. Estimaci´on de la Tendencia . . . . . . . . . . . . . . . . . . . . . . . . . .

2.3. Pron´osticos con base en la Tendencia . . . . . . . . . . . . . . . . . . . . .

2.4. Caso de Estudio: Pron´ostico de Ventas al Menudeo . . . . . . . . . . . . .

2.4.1. Descripci´on de los Datos . . . . . . . . . . . . . . . . . . . . . . .

2.4.2. Programa en R . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2.5. Problemas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3. Recursos en R Para An´alisis de Tendencia

3.1. Introducci´on . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3.2. Metodologias de Descomposici´on con Base en Suavizadores . . . . . . . .

3.2.1. Regresi´on Local Loess . . . . . . . . . . . . . . . . . . . . . . . .

3.2.2. STL: M´etodo de Descomposici´on Basada en Regresi´on Loess . . .

3.3. Filtros Lineales, Medias M´oviles y Suavizadores . . . . . . . . . . . . . .

3.4. La Funci´on Filter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3.5. Ejemplo de An´alisis de Tendencia con Suavizadores . . . . . . . . . . . . .

3.6. Notas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4. Modelado y Pron´ostico incluyendo la Componente Estacional

4.1. Definiciones . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4.2. Estimaci´on de la Componente Estacional . . . . . . . . . . . . . . . . . . .

4.2.1. Modelos con Variables Indicadoras . . . . . . . . . . . . . . . . .

4.2.2. Modelos con Funciones Trigonom´eticas . . . . . . . . . . . . . . .

4.3. Procedimientos Recursivos de Estimaci´on para Diagnosticar Estabilidad Estructural en Modelos de Pron´osticos . . . . . . . . . . . . . . . . . . . . .

4.4. M´etodos de Filtrado para la Componente Estacional . . . . . . . . . . . . .

vii 4.4.1. Algoritmo de Brockwell y Davis . . . . . . . . . . . . . . . . . . .

4.4.2. M´etodo de Holt-Winters . . . . . . . . . . . . . . . . . . . . . . .

4.5. Notas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4.6. Problemas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5. Validaci´on de los Supuestos sobre los Errores

5.1. Series Estacionarias en Covarianza . . . . . . . . . . . . . . . . . . . . . .

5.1.1. Estimaci´on de las funciones de Autocovarianza y Autocorrelaci´on. .

5.2. Pruebas de Incorrelaci´on . . . . . . . . . . . . . . . . . . . . . . . . . . .

5.2.1. Prueba Ljung-Box (LB) . . . . . . . . . . . . . . . . . . . . . . .

5.2.2. Prueba Durbin-Watson (DW) . . . . . . . . . . . . . . . . . . . . .

5.3. Alternativas cuando los Residuos Estructurales muestran Autocorrelaci´on .

6. Modelos ARMA para la Componente Aleatoria

6.1. Introducci´on . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6.2. Procesos de Medias M´oviles de orden q, MA(q) . . . . . . . . . . . . . . .

6.2.1. Propiedades . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6.2.2. Condici´on de Invertibilidad del Proceso MA(q) . . . . . . . . . . .

6.2.3. Funci´on fac parcial de un Proceso MA(q) invertible . . . . . . . . .

6.2.4. Implementaci´on en R . . . . . . . . . . . . . . . . . . . . . . . . .

6.3. Procesos Autorregresivos de Orden p, AR(p) . . . . . . . . . . . . . . . . . 100 6.3.1. Condici´on Suficiente para que un AR(p) sea Estacionario en Covarianza100 6.3.2. Algunas Propiedades de los Procesos AR(p) Estacionarios . . . . . 101 6.3.3. Proceso AR(1) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 6.4. Procesos Autoregresivos y de Medias M´oviles ARMA(p,q) . . . . . . . . . 107 6.4.1. Propiedades de los Modelos ARMA . . . . . . . . . . . . . . . . . 109

viii 6.4.2. Librer´ıas para identificaci´on, estimaci´on y pron´osticos de modelos ARMA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111 6.4.3. Identificaci´on de modelos ARMA . . . . . . . . . . . . . . . . . . 111 6.4.4. Estimaci´on de modelos ARMA . . . . . . . . . . . . . . . . . . . 112 6.5. Modelos ARMA Estacionales, SARMA . . . . . . . . . . . . . . . . . . . 114 7. Ra´ıces Unitarias y Tendencias Estoc´asticas (ARIMA)

119

7.1. Introducci´on . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119 7.2. MOdelos ARIMA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119 7.3. Ra´ıces Unitarias, Estimaci´on y Prueba de Hip´otesis . . . . . . . . . . . . . 127 7.3.1. Prueba Dickey-Fuller (DF) . . . . . . . . . . . . . . . . . . . . . . 129 7.3.2. Prueba Dickey-Fuller Aumentada . . . . . . . . . . . . . . . . . . 130 8. Ra´ıces Unitarias Estacionales y Estacionalidad Estoc´astica (SARIMA)

137

8.1. Modelos SARIMA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137 8.1.1. Modelo SARIMA

. . . . . . . . . . . . . . . . . . . . . . . . . . 137

8.2. Pruebas de Ra´ız Unitaria Estacional . . . . . . . . . . . . . . . . . . . . . 142 8.2.1. La prueba HEGY . . . . . . . . . . . . . . . . . . . . . . . . . . . 143 8.2.2. La prueba Canova-Hansen . . . . . . . . . . . . . . . . . . . . . . 150 A. Datos de Series de Tiempo

155

A.1. Series con Tendencia y Estacionalidad . . . . . . . . . . . . . . . . . . . . 155 A.1.1. Ejemplo 1. Cinco Series . . . . . . . . . . . . . . . . . . . . . . . 155 A.1.2. Ejemplo 2. Consumo domiciliario agua potable en Medell´ın . . . . 157

´Indice de figuras

2.1. Descomposici´on de la Serie Yt . . . . . . . . . . . . . . . . . . . . . . . .

2.2. Ajuste Log´ıstico a Poblaci´on de Medell´ın . . . . . . . . . . . . . . . . . .

2.3. Pron´osticos (—-) versus Observados (-o-o-) . . . . . . . . . . . . . . . . .

3.1. Suavizamiento Loess de N´umero de Licencias de Construcci´on en Medell´ın, 1986-2003 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3.2. Suavizamiento STL de N´umero de Licencias de Construcci´on en Medell´ın, 1986-2003 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3.3. Precios de la Acci´on de Microsoft Suavizados con SES . . . . . . . . . . .

3.4. Suavizamientos SES, Holt, Media Movil unilateral 30 del IGBC, 2009-2010

4.1. Average Monthly Temperatures at Nottingham, 1920-1939 . . . . . . . . .

4.2. Ajuste con Variables Indicadoras . . . . . . . . . . . . . . . . . . . . . . .

4.3. Estimadores Recursivos para el Modelo de Producci´on de Cemento

. . . .

4.4. Regi´on de Rechazo de la Prueba CU SU M . . . . . . . . . . . . . . . . . .

4.5. Resultados prueba CUSUM gr´afica . . . . . . . . . . . . . . . . . . . . . .

4.6. Punto de Quiebre 1974-Q4 en la Serie y en los Residuos . . . . . . . . . .

x 4.7. Ajuste con el M´etodo de Holt-Winters . . . . . . . . . . . . . . . . . . . .

4.8. Pron´osticos Varios M´etodos . . . . . . . . . . . . . . . . . . . . . . . . .

5.1. Ejemplo Funci´on de Autocovarianza . . . . . . . . . . . . . . . . . . . . .

5.2. Ejemplo Funci´on de Autocorrelaci´on. . . . . . . . . . . . . . . . . . . . .

5.4. Fac Te´orica de un Ruido Blanco. . . . . . . . . . . . . . . . . . . . . . . .

5.3. Ejemplos de Funciones de Autocorrelaci´on de un Proceso Estacionario . . .

5.5. Ejemplos de Fac y Bandas de Bartlett . . . . . . . . . . . . . . . . . . . .

5.6. Fac de los Residuos Estructurales Serie Cementos . . . . . . . . . . . . . .

5.7. Densidad χ2m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5.8. Fac Muestral de un AR(1). . . . . . . . . . . . . . . . . . . . . . . . . . .

5.9. Residuales estructurales versus Y estimada . . . . . . . . . . . . . . . . .

6.1. Funci´on de Autocorrelaci´on. . . . . . . . . . . . . . . . . . . . . . . . . .

6.2. FAC muestral de un M A(3). . . . . . . . . . . . . . . . . . . . . . . . . .

6.3. F AC y F ACP de un M A(3). . . . . . . . . . . . . . . . . . . . . . . . .

6.4. F AC y F ACP del Ejemplo. . . . . . . . . . . . . . . . . . . . . . . . . . 100 6.5. C´ırculo Unitario

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101

6.6. FAC de Yt = ϕYt−1 + εt . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 6.7. FACP Muestral de AR(p). . . . . . . . . . . . . . . . . . . . . . . . . . . 105 6.8. Fac Muestral de ARM A(p, q). . . . . . . . . . . . . . . . . . . . . . . . . 112 6.9. Posibles submodelos (1,1)=AR, (1,2),(2,1)=ARMA, (2,2)=SARMA . . . . 116 7.1. Trayectorias de Marchas Aleatorias . . . . . . . . . . . . . . . . . . . . . 121 7.2. Pron´osticos de la Serie PNB-USA. . . . . . . . . . . . . . . . . . . . . . . 124 7.3. (a) : Serie log(Precio) del Yen en USD, (b) Su fac, (c): La serie diferenciada, (d),(e): fac y facp de la serie diferenciada . . . . . . . . . . . . . . . . . . 124

xi 7.4. Pron´osticos usd/yen con ARIMA(3,1,2)(continua) y Tendencia Lineal+AR(2) (punteada) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125 7.5. Serie USD por Libra . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126 7.6. FAC y FACP Muestral de un ARIMA(p,1,q).

. . . . . . . . . . . . . . . . 128

8.1. Trayectoria de un modelo SARIMA(1,1,0)(1,1,0)[12] . . . . . . . . . . . . 139 8.2. Trayectoria de un modelo SARIMA(0,1,1)(0,1,1)[12] . . . . . . . . . . . . 141 8.3. Serie de Producci´on Trimestral de Cemento, Australia . . . . . . . . . . . . 146 8.4. Serie Diferenciada de Producci´on Trimestral de Cemento . . . . . . . . . . 148 8.5. Pron´osticos a 8 trimestres de la Producci´on Trimestral de Cemento . . . . . 149 8.6. Comparaci´on de los Ajustes y Pron´osticos. Panel superior: Modelo Descomposicion + AR, Panel inferior: SARIMA . . . . . . . . . . . . . . . . . . . 150

xii

CAP´ITULO

Introducion ´ al Lenguaje R

1.1. Introducci´on al Curso El curso de Series de Tiempo consiste de dos temas: el modelo de descomposici´on y el modelo ARIMA-SARIMA. El objetivo es desarrollar modelos de series de tiempo que permitan el c´alculo de pron´osticos. El curso est´a organizado con base en el texto de Diebold [1999]. Existen otros textos que desarrollan la teor´ıa de series de tiempo utilizando el R como el software b´asico, por ejemplo, Cowpertwait and Metcalfe [2009], Shumway and Stoffer [2005], Cryer and Chan [2008], Aragon [2008], Hyndman et al. [2008].

1.1.1.

Notas sobre Pron´osticos

Los pron´osticos se utilizan en forma constante en diversos campos: econom´ıa, finanzas, mercadeo, medio ambiente, ingenier´ıa. Una finalidad es proveer una gu´ıa para las decisiones que deben tomarse. Algunos ejemplos en estos campos pueden ser los siguientes. 1. Planeaci´on y Control de Operaciones. Por ejemplo, las decisiones de producci´on de un art´ıculo con base en los pron´osticos de ventas. Es posible por ejemplo, detectar una disminuci´on en la tendencia de ventas que conlleve a reducir la producci´on, o´ al contrario. 1

2 2. En Mercadeo la decisi´on de invertir en publicidad puede depender de pron´osticos de ventas. 3. En Econom´ıa. Las decisiones del Banco de la Rep´ublica, por ejemplo para el control de la inflaci´on, requieren el pron´ostico y el examen del comportamiento de ciertas variables macroecon´omicas, como el PIB, la tasa de desempleo, el IPC, las tasas de inter´es a distintos plazos, activas y pasivas. 4. En Econom´ıa los pron´osticos sobre ciclos econ´omicos de recesi´on y expansi´on para gu´ıa de aumento o´ disminuci´on de las tasas interbancarias. 5. En Planeaci´on. Los pron´osticos de demanda de viviendas, en los diferentes rangos, para orientar pol´ıticas de uso del suelo. 6. En Planeaci´on. Los pron´osticos de gastos sirven para definir los presupuestos futuros. 7. En Planeaci´on. Los pron´osticos de evoluci´on de casos de enfermedades o´ contingencias de salud. Aunque en casos como epidemias los modelos son series no lineales. 8. En Planeaci´on. El pron´ostico de consumo de energ´ıa el´ectrica domiciliaria es fundamental para las decisiones de generaci´on. 9. En Turismo. El pron´osticos de n´umero de turistas mensuales para determinar la demanda hotelera. 10. En Epidemiolog´ıa y Medio ambiente. La vigilancia de los niveles de contaminantes en el aire tiene como herramienta fundamental las series de tiempo. Pero adicionalmente el efecto de estos niveles sobre la salud.

1.2.

Introducci´on al R

1. Qu´e es R? http://www.r-project.org/. Es un software libre, de c´odigo abierto, para programaci´on orientada a objetos, dedicado a c´omputos estad´ısticos y financieros. 2. Qu´e es R-metrics? https://www.rmetrics.org/. Es un sistema de librer´ıas (toolboxs) de R para realizar finanzas computacionales e ingenier´ıa financiera. R posee aproximadamente 1600 librer´ıas sobre diversas aplicaciones. R-metrics posee 10 librer´ıas esplecializadas para manejo de riesgo, portafolios, valoraci´on, econometr´ıa,... 3. C´omo se obtiene R y se instala?.

3 a) Entrar a http://www.cran.r-project.org/. b) Seleccionar: “Download R for Windows” → “base”. c) Descargar el archivo .exe m´as reciente: (Download R 2.13.1 for Windows). d) Ejecutar el archivo .exe desde el escritorio. 4. Detalles de instalaci´on: R instala unas librer´ıas por defecto. El resto de las 1600 librer´ıas se puede instalar de manera inmediata dando click en (Paquetes→Instalar paquete(s)) donde se establece una conexi´on a Internet con los servidores de R en el mundo. A continuaci´on se escoge uno, por ejemplo, Colombia y luego aparece la lista de paquetes que se pueden instalar. Esta es la v´ıa para las librer´ıas de R-metrics. 5. Algunas librer´ıas: actuar: c´alculos actuariales. bbmle: estimaci´onpor m´axima verosimilitud. fBonds: c´alculos con bonos. fuzzyOP: l´ogica difusa. nnet: redes neuronales. 6. Algunos manuales: a) Entrar a http://www.cran.r-project.org/. b) Seleccionar: “Documentation”. c) Seleccionar: “Manuals” d) Seleccionar: “contributed documentation”, “Non-English Documents”. e) “R para Principiantes”, the Spanish version of “R for Beginners”, translated by Jorge A. Ahumada (PDF). f ) A Spanish translation of “An Introduction to R” by Andr´es Gonz´alez and Silvia Gonz´alez (PDF, Texinfo sources). g) “Metodos Estadisticos con R y R Commander” by Antonio Jose Saez Castillo (PDF, ZIP, 2010-07-08).

1.3.

Uso de R de Manera Interactiva versus Ejecuci´on desde un Programa

1.3.1.

Interactiva

Al ejecutar R aparece una ventana denominada R Console. Ventana interactiva. El s´ımbolo “>” indica que se espera un comando. > a = 7.5

#equivale a los dos siguientes

a a Crea un objeto num´erico “a” con el valor 7.5. > b = c(3,2,5) Crea un vector columna. > tb = t(b) Crea el vector traspuesto de “b”, es decir, como vector fila. > b1 = c(b,a) Concatena y crea un vector columna: b1 = (3, 2, 5, 7.5)0. > b2 = seq(1,length(b1),1) > b3 = cbind(b1,b2) Produce una matriz 4 × 2. > b4 = rbind(b1,b2) Produce una matriz 2 × 4. > (b34 = b3%*%b4) Produce producto matricial > b34 [1,] [2,] [3,] [4,]

[,1] [,2] [,3] [,4] 10.0 8 18.0 26.50 8.0 8 16.0 23.00 18.0 16 34.0 49.50 26.5 23 49.5 72.25

5 > h = matrix(0,4,4) Crea una matriz 4 × 4 de ceros. diag(h) = 1 Asigna unos a la diagonal de la matriz h, es decir, la vuelve una matriz identidad de 4 × 4. b34i = solve(b34 + h) Encuentra (b34 + h)−1 . bb = b34*h Multiplica elemento a elemento. > bb [1,] [2,] [3,] [4,]

1.3.2.

[,1] [,2] [,3] [,4] 10 0 0 0.00 0 8 0 0.00 0 0 34 0.00 0 0 0 72.25

Ejecuci´on desde un Programa

En lugar de o´ al mismo tiempo del manejo interactivo se pueden usar programas. Los programas en R tienen extenci´on .R o´ .r, por ejemplo, calculo.2.r. Los pasos para generar un programa se pueden describir como sigue. 1. Escoger un directorio base: Escogemos en la barra de men´u (Archivo→Cambiar dir...), luego se escoge o crea la ruta desde donde se quiera trabajar. Luego (Archivo→Nuevo script) y aparece otra ventana “Sin nombre - Editor R”. 2. En la barra men´u seleccionamos (Ventanas→Divididas verticalmente) para divir la pantalla. 3. Escribimos el programa. Por ejemplo “ejemplo1.r” # Ejemplo del uso de programas. # Muestra el uso de for, if, ifelse. b = 3 x = c(2,3,-4.2,5,10,7,3,2,3,0)

6 y = numeric(length(x)) for(i in 1:length(x)){ if(x[i] == b) y[i] = 0 else y[i] = 1 } z = ifelse(x==b,0,1) (cbind(y,z)) 4. Al terminar el c´odigo, en la pesta˜na (Archivo→Guardar como...), escribimos: ejemplo1.r 5. Para ejecutar el programa en la pesta˜na (Editar→Ejecutar todo) Suponga que se inicia una sesi´on nueva en R y se corre en la consola > source("ejemplo1.r") autom´aticamente se ejecutan los comandos en el programa. En la memoria quedan las variables generadas: b, x, y, z. > ls.str() muestra lo que hay en la memoria. > z = NULL elimina la variable z. Una parte importante de la programaci´on con R es escribir funciones. Para hacerlo se escribe el programa en el editor. Por ejemplo. # Funci´ on para calcular el estad´ ıstico t para la igualdad de medias prueba.t2 = function(y1,y2){ n1 = length(y1) n2 = length(y2) yb1 = mean(y1) yb2 = mean(y2) s1 = var(y1) s2 = var(2) s = ((n1-1)*s1 + (n2-2)*s2)/(n1 + n2 -2) tst = (yb1 - yb2)/sqrt(s*(1/n1 + 1/n2)) return(tst) } Se salva el archivo con el mismo nombre de la funci´on, es decir, prueba.t2.r Para cargarlo en una sesi´on interactiva o´ en otro programa se corre: source("prueba.t2.r")

7 Ejemplo 1.3.1. # Prueba de igualdad de medias asumiendo varianzas iguales. source("prueba.t2.r") s1 = 2; s2 = 2 u1 = 4; u2 = 6 n1 = 100 n2 = 30 x1 = rnorm(n=n1,mean=u1,sd=s1) x2 = rnorm(n2,u2,s2) tst = prueba.t2(x1,x2) # C´ alculo del valor cr´ ıtico de una t con v=n1+n2-2 # Para alpha = 0.05/2 = 0.025 v = n1 + n2 - 2 tcrit = abs(qt(0.025,v)) abs(tst) < tcrit # Si TRUE no se rechaza H0:u1=u2

Concepto de data.frame (Hojas de Datos) Un objeto data.frame es una clase de objeto en R que es una colecci´on de columnas que contienen datos, que no tienen que ser del mismo tipo, con el mismo n´umero de filas. Las columnas pueden tener nombres. (ver Cap 6 de Manual o´ secci´on 2.2.2 de Farnsworth). Ejemplo 1.3.2. Si D se ha definido como data.frame, y contiene las columnas edad, genero, salario, entoces al ingresar ( > D\$edad ) en la consola, escribe los datos de la columna edad.

Concepto de list Una lista (list) es otro tipo de objeto en R, m´as general que un data.frame. En una lista se coloca un grupo de otros objetos que pueden ser: escalares, data.frames, cadenas de s´ımbolos. Las funciones en R generalmente retornan un objeto que puede ser un escalar,o tambi´en vector, o matr´ız, o un data.frame, o un list. Ejemplo 1.3.3. nombre = c("Cecilia", "Olga", "Beatriz") ab = 1:3

8 cd = 3:5 # crear un data.frame (B = data.frame(nombre=nombre, ab=ab, cd=cd)) B$nombre[2] # Olga # crear un list con un data.frame inclu´ ıdo R = list(B=B, ab=ab, cd=cd) R$B$nombre[1] # Cecilia

1.3.3.

Lectura de Datos

R posee varias funciones para lectura de datos. Incluso, tiene capacidad para leer archivos que se encuentran en un sitio Web determinado, dado que exista una conexi´on de internet.

La funci´on read.table La funci´on G = read.table("nombre") se utiliza para leer datos en R de tal forma que retorna un objeto data.frame, en este ejemplo indicado por la letra G, y donde “nombre” es el nombre del archivo en donde est´an los datos, por ejemplo, nombre = base.dat. Las extensiones pueden ser “.dat”, “.txt”, “.prn”. En particular, “.prn” es la adecuada cuando se quieren leer datos de excel. En este caso se guarda el archivo excel con el formato “.prn”, que corresponde a columnas separadas por espacios. No es conveniente salvar un archivo excel con extensi´on “.txt” porque este formato guarda controles de tabuladores invisibles que impiden la lectura de R. Si cada columna en el archivo de lectura de datos tiene un nombre entonces el comando se modifica colocando G = read.table("nombre", header = TRUE). Entonces conviene a˜nadir la instrucci´on attach(G), la cual hace que los nombres de las columnas pasen a ser los nombres de las variables activas en R. Ocasionalmente una de las columnas es la fecha de la observaci´on, con el formato d´ıa-mes-a˜no, por ejemplo, “23/10/2010”. En este caso es conveniente colocar la instrucci´on en la forma G = read.table("nombre", header = TRUE, stringsAsFactors=FALSE) El efecto es que toma los datos de las fechas como alfanum´ericos y no como factores, que es un tipo de datos para utilizar en dise˜no de experimentos. Ejemplo 1.3.4. En este ejemplo se muestra c´omo utilizar la funci´on read.table para leer datos de un archivo en una p´agina web. El archivo tiene dos columnas con los nombres: fecha y x. Los datos de la variable fecha son alfanum´ericos.

9 #---------------------------------------archivo = "http://www.unalmed.edu.co/˜ndgirald/ Datos/Datos%20curso%20Series%20II/fechaydatos.prn" G = read.table(archivo, header = TRUE, stringsAsFactors=FALSE) attach(G) np = length(x) fechas = strptime(fecha, "%d/%m/%y") plot(fechas,x, xaxt="n", panel.first = grid(),type=’l’,ylab=’’) axis.POSIXct(1, at=seq(as.Date(fechas[1]),as.Date(fechas[np]), "months"), format="%m/%y") axis.POSIXct(1, at=seq(as.Date(fechas[1]),as.Date(fechas[np]), "years"), labels = FALSE, tcl = -0.2)

Otras funciones de lectura y escritura Hay otras funciones de lectura, como G = read.csv("nombre", header = TRUE), que lee archivos excel “delimitados por comas”, de extensi´on .csv. Y funciones para escribir datos como la funci´on G = write.matrix(data.frame, archivo). #---------------------------------------D = data.frame(fechas=fecha,x=x) require(MASS) write.matrix(D,"h.dat")

1.4. Repaso del Modelo de Regresi´on Lineal El modelo de regresi´on lineal m´ultiple es b´asico para algunos de los modelos de series de tiempo. Suponga que se tienen unas muestras de las variables aleatorias Y, X1 , X2, Yt , X1,t, X2,t,

t = 1, . . . , T.

Adem´as suponga que εt ∼ iid N (0, σ2) es una muestra de una Normal. Si se cumple que Yt = β0 + β1 X1,t + β2 X2,t + εt ,

(1.1)

entonces se dice que las variablesY, X1, X2 satisfacen un modelo de regresi´on lineal. Donde Yt es la variable dependiente y X1,t, X2,t son las variables explicativas y el t´ermino εt es

10 el error o´ residuo. Se asume que las variables X1,t, X2,t est´an dadas. Los supuestos del modelo lineal son 1. Las εt no est´an correlacionadas. 2. Las εt tienen varianzas constante. 3. εt es independiente de las Xi,t. Lo cual implica a) E(εt |X1,t, X2,t) = E(εt ) = 0, luego E(Yt|X1,t, X2,t) = β0 + β1 X1,t + β2 X2,t. b) V ar(Yt |X1,t, X2,t) = V ar(εt|X1,t, X2,t) = V ar(εt) = σ 2 . donde σ 2 se denomina la varianza del error, y es constante. 4. El rango de la matriz X es 3, es decir, no existen constantes a, b tal que aX1,t +bX2,t = 1, es decir, las columnas de X son linealmente independientes. El modelo se dice que es lineal porque Yt se relaciona linealmente con X1,t, X2,t, es decir, Yt es una combinaci´on lineal de X1,t, X2,t. Una relaci´on, por ejemplo Yt =

1 1

+ eβ0 +β1 X1,t +β2 X2,t

+ εt ,

(1.2)

es no lineal.

1.4.1.

Estimaci´on de los Par´ametros

Los par´ametros de los modelos (1.1) y (1.2), son: β0 , β1 , β2 , σ 2 y se busca estimarlos a partir de la muestra de Yt , X1,t, X2,t. Algunos m´etodos de estimaci´on son: 1. M´ınimos cuadrados ordinarios. (MCO) 2. M´ınimos cuadrados robustos. (MCR) 3. M´axima verosimilitud. (MLE) 4. M´ınimos cuadrados no lineales. 5. Momentos.

11 En el caso MCO, los estimadores de los coeficientes βi se definen como el vector βˆ = (βˆ0 , βˆ1, βˆ2 ) que minimiza la distancia G(β) =

T X t=1

2 Yt − (β0 + β1 X1,t + β2 X2,t) ,

(1.3)

ˆ ≤ G(β). Y se puede escribir es decir, βˆ cumple que ∀β, G(β) ˆ = argmin G(β) β β

Minimizar la funci´on objetivo G(β) es resolver el problema de estimaci´on por MCO. El ˆ se calcula mediante una f´ormula. Si se define la matriz de dise˜no XT ×3 como estimador β 

  X=  

1 X1,1 1 X1,2 .. .. . . 1 X1,T

X2,1 X2,2 .. . X2,T





    , Y =     

Y1 Y2 .. . YT





    , ε =     

ε1 ε2 .. . εT



  ,  

entonces Y = Xβ + ε es el modelo en forma matricial y se puede comprobar que βˆ = (X 0X)−1 X 0 Y es el estimador MCO. ˆ = β, de m´ınima varianza, Este estimador βˆ es un estimador lineal insesgado, es decir, E(β) ? ? osea, V ar(βˆi) ≤ V ar(βˆi ) para βˆi cualquier estimador diferente del de MCO. El problema de estimaci´on por MCO no requiere que εt ∼ iid N (0, σ2). En cambio, la estimaci´on MLE s´ı. El supuesto εt ∼ iid N (0, σ 2) implica que los βbi ∼ N (·, ·) y se pueden definir los estad´ısticos de ajuste como la t y la F. Note que una vez definidos los βˆi, se define los residuos estimados como εˆt = Yt − βˆ0 + βˆ1 X1,t + βˆ2 X2,t ,

t = 1, . . . , T.

(1.4)

Y se definenen los pron´osticos dentro de la muestra (interpolar) como Yˆt = βˆ0 + βˆ1 X1,t + βˆ2 X2,t,

t = 1, . . . , T,

(1.5)

y fuera de la muestra (extrapolar) como YˆT +h = βˆ0 + βˆ1 X1,T +h + βˆ2 X2,T +h , donde X1,T +h , X2,T +h hay que proveerlos.

h = 1, . . . , m,

(1.6)

1.4.2.

Pruebas de Ajuste del Modelo

Pruebas de Significaci´on de los par´ametros βi Si εt ∼ iidN (0, σ2) se cumple que las βˆi ∼ N (βi, σ 2ˆ ), donde σ ˆβˆi = Sβˆi = error est´andar. βi Adem´as se cumple que t = βˆi /S ˆ se distribuye tT −k , donde k es el n´umero de par´ametros βi

βi , y es el estad´ıstico de la prueba de significancia de los par´ametros βi . Esta prueba se describe as´ı. Para i = 1, . . ., k, 1. H0 : βi = 0 versus H1 : βi 6= 0. 2. Estad´ıstico de prueba ti = βˆi /Sβˆi . 3. Si H0 es cierta entonces ti ∼ tT −k . 4. Decisi´on. Se plantean 3 formas equivalentes. 1 Si |ti | > 1.96 entonces se rechaza H0 : bi = 0 con un nivel de significancia de 5 %, es decir, hay una probabilidad de 5 % de rechazar H0 siendo cierta. 2 Si 0 6∈ [βˆi − 1.96Sβˆi , βˆi + 1.96Sβˆi ] se rechaza a un nivel de significancia de 5 %. 3 Defina V alorp = P (|tT −k | > |tobs |) = P (tT −k > |tobs |)+P (tT −k < −|tobs |). Si V alorp < 0.05 se rechaza H0 a un nivel de significancia del 5 %, de lo contrario no se rechaza H0 .

Prueba F Se define la suma de cuadrados total SST = SSR + SSE como la suma de cuadrados regresi´on m´as suma de cuadrados de errores, dada por SST =

T X (Yt − Y¯ )2

(1.7)

t=1

T T X X 2 ˆ ¯ (Yt − Yˆt )2 (Yt − Y ) + = t=1

(1.8)

t=1

= SSR + SSE.

Se define el estad´ıstico F para probar la hip´otesis H0 : β1 = β2 = · · · = βk−1 = 0 versus H1 : no(H0 )

(1.9)

13 donde el estad´ıstico est´a dado por SSR H F = k − 1 ∼0 Fk−1,T −k . SSE T −k

(1.10)

Si se rechaza H0 , las variables X1,t, . . . , Xk−1,t tienen capacidad para predecir los valores de Yt . Si F > Fk−1,T −k,α se rechaza H0 . Ejemplo 1.4.1. Con k = 3, T = 48, α = 0, 05, se tiene que Fk−1,T −k,α = F2,45,0.05 = 3.2, luego si F = 30.89. Como F > 3.2 se rechaza H0 , es decir, X1,t y X2,t tienen capacidad predictiva sobre Yt .

1.4.3.

Estad´ısticos de Ajuste y de Selecci´on de Modelos

Se define el estimador de la varianza del error como T T 1 X 2 1 X 2 ˆ (Yt − Yt ) = εˆt , σ ˆ = T −k T −k 2

t=1

(1.11)

t=1

donde k es el n´umero de β’s en el modelo. Se denomina “el error cuadr´atico medio estimado” o´ MSE estimado. Se define el estad´ıstico de ajuste R cuadrado, R2 , como 2

R = 1 − PT

ˆ2t t=1 ε

t=1 (Yt

− Y¯ )2

(1.12)

como el porcentaje de la varianza de Yt que es explicado (atribuible) por las variables explicativas. Observe que k no aparece. R2 ≥ 0.7 buena

R2 = 0.6 buena pero no mucho R2 = 0.4 regular-mala R2 = 0.2 desechar, seguramente la prueba F y las pruebas t no rechazan ¯ 2 como Se define el R cuadrado ajustado, R ¯ 2 = 1 − (T − 1)SSE . R (T − k)SST

(1.13)

14 ¯ 2 se puede escribir en terminos de R2 de la Note que R2 = 1 − SSE/SST , por tanto, R siguiente forma ¯ 2 = 1 − T − 1 (1 − R2 ). R (1.14) T −k Si k %, entonces T − k & y n´umero de par´ametros.

T −1 T −k

¯ 2 &. Es decir, R ¯ 2 penaliza aumentar el %, luego R

Se define el criterio de informaci´on de Akaike AIC PT ˆ2t (T − k)ˆ σ 2 2k/T 2k/T t=1 ε AIC = e = e . (1.15) T T P donde Tt=1 εˆ2t = (T − k)ˆ σ 2 . El AIC es un estimador de σ 2 pero penalizado por el n´umero de grados de libertad, es decir, aumenta cuando k aumenta. El AIC tambi´en se define como el logaritmo de (1.15). De la definici´on del criterio AIC se observa: 2k

Si k %, entonces e T % luego AIC %. Si M SE &, entonces AIC &. ((Como en el caso de S 2 muchos de los criterios m´as importantes para la selecci´on de modelos de pron´ostico tienen la forma de factor de penalizaci´on multiplicado por M SE)) Diebold [1999, p´ag. 73] Se define el criterio de informaci´on de Schwarz BIC PT ˆ2t k/T t=1 ε BIC = T T

(1.16)

El BIC tambi´en se define como el logaritmo de (1.16). La regla para utilizar AIC y BIC es que para escoger entre varios modelos de regresi´on con respecto a la misma variable dependiente, todos lineales, o´ todos no lineales pero con la misma estructura, es decir, modelos anidados, se escoge el de menor AIC o´ menor BIC. Diebold [1999, p´ag. 74–75], introduce dos criterios para elegir entre AIC y BIC. 1. Consistencia. 2. Eficiencia asint´otica.

15 Definici´on 1.4.1 (Consistencia). Si se examinan varios modelos, y el modelo generador de los datos (MGD) est´a inclu´ıdo, el criterio se dice consistente si a medida que T aumenta, la probabilidad de que el criterio seleccione el MGD tiende a 1. Resultado: S 2 , AIC no son consistentes pero BIC s´ı es. Definici´on 1.4.2 (Eficiencia Asint´otica). Un criterio es eficiente asint´oticamente si a medida que T aumenta, elige varios modelos cuyas varianzas de error de pron´ostico a un paso tienden a la que se obtendr´ıa con el MGD. Donde el error de pron´ostico a un paso est´a dado por YT +1 − YˆT +1 Resultado: AIC es eficiente asint´oticamente pero BIC no.

((Muchos autores recomiendan usar el modelo m´as parsimonioso que selecciona el BIC en igualdad de circunstancias)) Diebold [1999, p´ag. 75].

1.4.4.

M´ınimos Cuadrados Nolineales

En algunos casos es necesario ajustar un modelo no lineal. R posee funciones que permiten ajustar por m´ınimos cuadrados no lineales modelos de la forma Y = g(X; β) +

(1.17)

N´otese que debe ser un modelo con errores aditivos. La funci´on en R para m´ınimos cuadrados no lineales es nls() y se aplica de la misma forma que lm(). Por ejemplo, para el modelo no lineal (1.2) 1 + εt , (1.18) Yt = β +β X1,t +β2 X2,t 0 1 1+e Se puede estimar con la funci´on nls(), con la instrucci´on m = nls(Y ˜ 1/(1+exp(w+a*X1+b*X2)), start=list(w=0.1, a=1, b=1)) N´otese que la funci´on requiere unos valores iniciales para los par´ametros. Con la transformaci´on Z = log(1/Y − 1) se puede linearizar el modelo y colocar Z = β0 + β1 X1,t + β2 X2,t. Con este modelo se podr´ıan obtener estimadores iniciales de los par´ametros. En el cap´ıtulo siguiente hay un ejemplo concreto de estimaci´on de este modelo.

1.5.

Ejemplo de Regresi´on con R

En este ejemplo se muestran algunas de las capacidades de R. En particular: 1) La lectura de datos en un sitio de la web a partir de la direcci´on url. 2) Los estad´ısticos para escogencia del mejor modelo. 3) Gr´aficas de los modelos ajustados. Los datos son de tres variables Y, X1, X2, para un modelo de regresi´on (1.1). Utilizamos datos de un ejemplo de la p´agina web http://www.stat.umn.edu/geyer/5102/examp/reg.html Se examinan tres modelos lineales: Yt = β0 + β1 X1,t + εt ,

(1.19)

Yt = β0 + β1 X1,t + β2 X2,t + εt ,

(1.20)

Yt = β1 X1,t + β2 X2,t + εt ,

(1.21)

Los resultados aparecen en el Cuadro 1.1

#-------------------------------------------------------------# programa en R library(xtable) # para reportar tablas en LaTex para informes library(fRegression) # para estadisticos especiales de regresion archivo = "http://www.stat.umn.edu/geyer/5102/data/ex5-3.txt" D = read.table(archivo,header=T) attach(D) m1 = lm(y ˜ x1) summary(m1) print(xtable(m1),digits=4) m2 = lm(y ˜ x1 + x2) summary(m2) print(xtable(m2),digits=4) m3 = lm(y ˜ -1 + x1 + x2) summary(m3) print(xtable(m3),digits=4)

Tabla 1.1: Resultados de los Modelo Regresi´on

(Intercept) x1 (Intercept) x1 x2 x1 x2

Estimate 2.3778 1.6979 Estimate 0.4342 1.4179 0.6743 Estimate 1.4623 0.7174

Std. Error 0.8797 0.1798 Std. Error 0.9094 0.1719 0.1688 Std. Error 0.1435 0.1415

t value 2.70 9.45 t value 0.48 8.25 3.99 t value 10.19 5.07

Pr(>|t|) 0.0095 0.0000 Pr(>|t|) 0.6352 0.0000 0.0002 Pr(>|t|) 0.0000 0.0000

Las Tablas Anova para los modelos m2, m3 son las siguientes. #-----------------------------------------------------anova(m2) Residual standard error: 0.9722 on 47 degrees of freedom Multiple R-squared: 0.7388, Adjusted R-squared: 0.7277 F-statistic: 66.48 on 2 and 47 DF, p-value: 1.987e-14 #-----------------------------------------------------anova(m3) Residual standard error: 0.9644 on 48 degrees of freedom Multiple R-squared: 0.9922, Adjusted R-squared: 0.9919 F-statistic: 3061 on 2 and 48 DF, p-value: < 2.2e-16 Y los valores de AIC aparecen a continuaci´on. #-------------------------(c(AIC(m1),AIC(m2),AIC(m3))) 156.5899 143.9822 142.2242

18 Se podr´ıa escoger el modelo m2 ya que los valores del R cuadrado ajustado para el modelo 3 no son v´alidos por ser un modelo sin intercepto. Adem´as, los valores estimados de la variable Yt , dentro de la muestra, son muy similares para ambos, el modelo 2 y el modelo3.

1.6.

Manejo de Series de Tiempo con R

En esta secci´on consideramos una serie de tiempo Tt , t = 1, 2, . . ., T , como un vector columna. En R se ingresar´ıa inicialmente como un objeto tipo vector num´erico. Pero para los an´alisis en los cap´ıtulos siguientes se requiere redefinir la serie como un objeto de clase , mediante la funci´on ts().

CAP´ITULO

Modelado y Pronostico ´ de la Tendencia

Una serie de tiempo es una sucesi´on de variables aleatorias ordenadas de acuerdo a una unidad de tiempo, Y1 , . . . , YT . Por ejemplo, la concentraci´on en el aire de di´oxido de azufre en ppm (partes por mill´on), (100ppm = 262mg/m3), medida semanalmente en un punto espec´ıfico, es importante para monitorear la calidad del aire en una ciudad. Pueden mencionarse varios motivos por los cuales es u´ til analizar las series de tiempo. Por ejemplo, para obtener pron´osticos de valores futuros de la serie, con el fin de ayudar a tomar decisiones que tienen consecuencias importantes. O para entender mejor el mecanismo de generaci´on de los datos, que puede no ser claro inicialmente en una investigaci´on.

2.1. El Modelo Aditivo de Componentes de Series de Tiempo Dada una serie Yt , t = 1, . . . , T , el Modelo Aditivo de Componentes consiste en asumir que Yt se puede descomponer en tres componentes: Yt = Tt + St + εt ,

(2.1)

donde Tt es la tendencia, St es la componente estacional y εt es la componente de errores. Las componentes Tt y St son funciones de t determin´ısticas. Su evoluci´on es perfectamente predecible. 19

20 La componente Tt en algunos casos tambi´en puede ser una componente estacional, pero de baja frecuencia, o, equivalentemente, una componente con per´ıodo muy grande. Por ejemplo, en una serie diaria, St puede tener per´ıodo 30 d´ıas, y Tt per´ıodo 360 d´ıas. El Modelo Multiplicativo consiste en asumir que Yt se puede descomponer en tres componentes: Yt = TtSt exp εt . (2.2)

5 4

6 4

En la Figura (2.1) siguiente se muestra la idea de la descomposici´on. Al superponer las series en los paneles (a), (b) y (c) se obtiene la serie en el panel (d).

60 Time

(a) Tendencia

100

Time

(b) Patr´on Estacional

100

Time

100

Time

(d) Serie

Figura 2.1: Descomposici´on de la Serie Yt El an´alisis consiste en modelar y estimar Tt y St y luego extraerlas de Yt para obtener εˆt = Yt − Tˆt − Sˆt . La serie εˆt se modela y estima para finalmente reconstru´ır Yt , Yˆt = Tˆt + Sˆt + εˆt , y poder realizar el pron´ostico YˆT +h = TˆT +h + SˆT +h + εˆT +h , utilizando la informaci´on disponible Y1 , . . ., YT con h = 1, 2, . . ., m. Sin embargo, puede suceder que la serie εˆt sea incorrelacionada, es decir, Corr(ˆ εt , εˆt+s ) = 0, para s 6= 0. En este caso εˆT +h = 0, ∀h > 0. Definici´on 2.1.1 (Tendencia). Se define como una funci´on Tt de t que describe la evoluci´on lenta y a largo plazo del nivel medio de la serie. La funci´on Tt depende de par´ametros, que deben estimarse. Modelos: Una lista de posibles modelos para Tt es: Lineal Tt = β0 + β1 t,

(2.3)

21 Cuadr´atico Tt = β0 + β1 t + β2 t2 ,

(2.4)

Tt = β0 + β1 t + β2 t2 + β3 t3 ,

(2.5)

Tt = exp(β0 + β1 t).

(2.6)

C´ubico

Exponencial

Log´ıstico Tt =

β2 , 1 + β1 exp(−β0 t)

(2.7)

En la tendencia cuadr´atica se puede observar: 1. Si β1 , β2 > 0, Tt es mon´otona creciente. 2. Si β1 , β2 < 0, Tt es mon´otona decreciente. 3. Si β1 > 0 y β2 < 0, Tt es c´oncava hacia abajo (o c´oncava). 4. Si β1 < 0 y β2 > 0, Tt es c´oncava hacia arriba (o convexa). Definici´on 2.1.2. El modelo Logar´ıtmico Lineal o´ Log-Lineal se define como ln Yt = β0 + β1 t + εt .

(2.8)

Corresponde a un modelo con tendencia lineal para el logaritmo de Yt . En (2.8) al tomar exponencial se tiene Yt = exp(β0 + β1 t + εt ), que es similar al modelo con tendencia exponencial (2.6), Yt = exp(β0 + β1 t) + εt . Sin embargo, son modelos diferentes y se estiman por m´etodos diferentes. Ejemplo 2.1.1. En el Cuadro 2.1 se muestran los datos del total de habitantes en Medell´ın, seg´un los Censos desde 1905 hasta 2005, seg´un el DANE (1 ). En Poveda [1982] se plan1

http://es.wikipedia.org/wiki/Demograf´ıa de Medell´ın

Tabla 2.1: Poblaci´on de Medell´ın censos 1905 - 2005

1 2 3 4 5 6 7 8 9 10 11

A˜no 1905 1912 1918 1928 1938 1951 1964 1973 1985 1993 2005

Poblaci´on 59.810 70.550 79.150 120.040 168.270 358.190 772.890 1.077.250 1.468.090 1.630.010 2.223.080

te´o ajustar los datos de desde 1938 hasta 1973 mediante un modelo de componentes con tendencia log´ıstica, (2.7). β2 + εt . 1 + β1 exp(−β0 t)

(2.9)

1000 500 0

mil ones de hab

1500

2000

Yt =

1920

1940

1960

1980

Año

Figura 2.2: Ajuste Log´ıstico a Poblaci´on de Medell´ın

2000

2.2. Estimaci´on de la Tendencia En este cap´ıtulo se introduce la estimaci´on de la tendencia mediante modelos de regresi´on lineal y no lineal. Son modelos param´etricos. En el cap´ıtulo siguiente se introducen modelos no param´etricos para estimar la tendencia, como los suavizadores, los filtros lineales y no lineales y las medias m´oviles. Hay otros m´etodos que no se consideran en este curso, por ejemplo, wavelets. En ocasiones la expresi´on “suavizar una serie” es equivalente a “extracci´on de la tendencia de una serie”, y ambas equivalen a la estimaci´on de la tendencia. Para la estimaci´on de los par´ametros β = (β0 , β1 , β2 )0 en los modelos lineales (2.3),(2.4), (2.5) y (2.8) se utiliza el m´etodo de m´ınimos cuadrados ordinarios (MCO). En este m´etodo el vector de par´ametros estimados βˆ es el vector que produce el valor m´ınimo de la suma de 2 P errores cuadrados. Es decir βˆ es el valor en el cual G(β) = Tt=1 Yt − Tt (β) toma el valor m´ınimo. βˆ = argmin G(β). (2.10) β

Para los modelos (2.6) y (2.7) se usa el m´etodo de m´ınimos cuadrados no lineales, que 2 P tambi´en minimiza la suma de errores cuadrados G(β) = Tt=1 Yt − Tt(β) , pero Tt (β) es una funci´on no lineal de β. El modelo LogLineal (2.8) es equivalente, algebr´aicamente, a Yt = exp(β0 + β1 t + εt ).

(2.11)

Sin embargo, este u´ ltimo modelo es no lineal y no coincide con el modelo exponencial,(2.6), Yt = exp(β0 + β1 t) + εt . Es posible estimar por m´ınimos cuadrados ordinarios el modelo LogLineal y utilizar los par´ametros estimados βˆ0 , βˆ1 como valores iniciales en la estimaci´on del modelo exponencial por m´ınimos cuadrados no lineales. Pero los par´ametros estimados en ambos modelos no necesariamente coinciden. Aunque la serie tenga una componente estacional St , Yt = Tt + St + εt , solamente consideramos un modelo de regresi´on entre Yt y Tt , tal que Yt = Tt + ηt, donde ηt es el t´ermino de error, de forma que ηt = St + εt . Por ejemplo. 1. En el caso lineal β = (β0 , β1 )0 , con Tt = β0 + β1 t se ajusta el modelo de regresi´on lineal: Yt = β0 + β1 t + ηt. 2. En el caso cuadr´atico β = (β0 , β1, β2 )0 , con Tt = β0 + β1 t + β2 t2 se ajusta el modelo de regresi´on lineal Yt = β0 + β1 t + β2 t2 + ηt. N´otese que en este caso hay que definir variable explicativa adicional t2 .

2.3.

Pron´osticos con base en la Tendencia

Suponga la serie con tendencia Yt = Tt +ηt , t = 1, . . . , T , con (ηt) una sucesi´on iid(0, σ 2). Los pron´osticos de Yt en los tiempos T + 1, T + 2, . . . , T + h, h ≥ 1 se definen como Definici´on 2.3.1. YbT +j = TbT +j , j = 1, . . . , h.

(2.12)

donde Tbt es la funci´on estimada de la tendencia. Por ejemplo, en el modelo lineal Yt = β0 + β1 t + εt , al reemplazar t por T + h se obtiene YT +h = βˆ0 + βˆ1 (T + h) + εˆT +h . Pero el pron´ostico ηˆT +h , puede o no ser cero. Es cero si los residuos εˆt , t = 1, . . . , T son independientes. Para decidir esto se realiza una prueba de incorrelaci´on y una prueba de normalidad. Puede suceder que los residuos estimados εˆt sean autocorrelacionados. En este caso el pron´ostico es YbT +j = TbT +j + εˆT +j . La definici´on general de pron´ostico, para una serie Yt , t ∈ Z, con base en la informaci´on Y1 , . . . , YT es una esperanza condicional, como sigue. Definici´on 2.3.2. YbT +j = E(YT +j |Y1 , . . . , YT ), j = 1, . . ., h.

(2.13)

Otros tipos de pron´osticos son

El pron´ostico por intervalo se obtiene si se asume que εt ∼ iid N (0, σ 2). Entonces se cumple que YˆT +h ± 1.96ˆ σ es un intervalo de confianza del 95 %. Si se ignora los errores en βˆ0 y βˆ1 el IC es βˆ0 + βˆ1 (T + h) ± 1.96ˆ σ, donde σ ˆ2 = P T 1 2 M SE = T −2 t=1 et El pron´ostico de densidad es YT +h ∼ N (YˆT +h, σ ˆ 2 ).

2.4.

Caso de Estudio: Pron´ostico de Ventas al Menudeo

Este caso est´a analizado en Diebold [1999, secci´on 4.5]. El objetivo aqu´ı es repetir el an´alisis utilizando R. Los Modelos a utilizar son los modelos lineal, cuadr´atico, c´ubico y exponencial. Se aplicar´an los criterios AIC y BIC para escoger el m´as adecuado.

2.4.1.

Descripci´on de los Datos Ventas al menudeo en USD a precios de USD en 1999 Periodicidad: Mensual Fechas: 01/1955 - 12/1994 N´umero de observaciones: 468 Datos ajustados estacionalmente. Es decir, si se tiene Zt = Tt + St + εt , se dice que Yt est´a ajustada estacionalmente o desestacionalizada si Yt = Zt − St = Tt + εt . Es decir, es una serie en la que se elimin´o la componente estacional. Usar el per´ıodo 01/1955 - 12/1993 para estimar los modelos (per´ıodo de entrenamiento) y el periodo 01/1994 - 12/1994 para examinar la eficiencia de los pron´osticos fuera de la muestra. Esta t´ecnica para verificar la eficiencia del modelo estimado se llama “validaci´on cruzada”.

2.4.2.

Programa en R

Lectura de Datos, Graficaci´on con fechas D = read.table("ventas_al_menudeo.dat",header=T) attach(D) # utiliza el nombre de las columnas como variables # # # # #

no hay variable con fechas : mensual 01/1955 hasta 12/1994 para 469 obs, des-estacionalizadas RTTR es el volumen de ventas en grandes almacenes en usd de 1999 La variable RTRR del archivo tiene datos faltantes NA al inicio y al final

y = na.omit(RTRR)/10000 # Convertir los datos en un objeto tipo ts y = ts(y,frequency=12,start=c(1955,01)) # generar un vector de fechas, clase ’Date’ fechas = seq(as.Date("1955/1/1"), length.out = length(y), by =

"months")

26 # grafica con fechas ts.plot(y,main="Ventas al menudeo en UDS de 1999") # otros comandos para graficar con fechas con mas detalle: mes-a˜ no np = length(y) ejex.mes = seq(fechas[1],fechas[np], "months") ejex.a˜ no = seq(fechas[1],fechas[np],"years") plot(fechas,y, xaxt="n", panel.first = grid(),type=’l’, ylab=’ventas.mes’, lwd = 2) axis.Date(1, at=ejex.mes, format="%m/%y") axis.Date(1, at=ejex.a˜ no, labels = FALSE, tcl = -0.2)

Ajuste de los Modelos utilizando Regresi´on Lineal Ajustar cuatro modelos: lineal, cuadr´atico, c´ubico y exponencial, mediante regresi´on lineal, utilizando la funci´on de R, lm().

# Generar datos para validacion cruzada: dejar el ultimo a˜ no T = length(y) yi = y[1:(T-12)] yf = y[(T-12+1):T] # Ajustar 4 modelos: lineal, cuadr´ atico, c´ ubico, log-lin t = seq(1:(T-12)) t2 = tˆ2 t3 = tˆ3 lyi = log(yi) # estimacion por minimos cuadrados ordinarios mod.lin = lm(yi˜t) mod.cuad = lm(yi˜t+t2) mod.cub = lm(yi˜t+t2+t3) mod.llin = lm(lyi˜t) # auxiliar para el exponencial

summary(mod.lin) summary(mod.cuad) summary(mod.cub)

Ajuste del Modelo Exponencial Lineal # paso 1) estimar el modelo auxiliar log - linear mod.llin = lm(lyi˜t) # paso 2) guardar los parametros del log-lineal b0.est = mod.llin$coefficient[1] b1.est = mod.llin$coefficient[2] # paso 3) guardar los datos en un data.frame Ds = data.frame(yi,t) # paso 4) usar la funcion nls mod.exp = nls(yi˜exp(beta0+beta1*t), data=Ds, start=list(beta0=b0.est, beta1=b1.est)) # paso 5) resultados summary(mod.exp)

Ajuste de los Modelos N´otese de los resultados del cuadro 2.2 que el modelo c´ubico no ajusta ya que el coeficiente de t3 no d´a significativo al nivel de 5 %. Se procede a calcular los estad´ısticos AIC y BIC.

Calcular los Estad´ısticos de Selecci´on del Modelo Para calcular los estad´ısticos de Selecci´on: R2 ajustado, MSE, log(AIC), log(BIC) se usar´a la funci´on medidas(), que los calcula para cada modelo. medidas = function(m,y,k){ # y = serie, m = modelo, k = numero parametros

Tabla 2.2: Ajuste de los Modelos Lineal, Cuadr´atico, C´ubico y Exponencial

(Intercept) t (Intercept) t t2 (Intercept) t t2 t3 (Intercept) t

Estimate -17537.673441 340.738727 18886.034127 -111.729690 0.938731 19386.588312 -124.127295 1.002967 -0.000089 9.390e+00 5.769e-03

Std. Error 1503.384304 5.405165 383.923338 3.678499 0.007390 513.151618 9.210396 0.044379 0.000061 1.432e-02 3.536e-05

t value -11.665463 63.039465 49.192201 -30.373717 127.020988 37.779455 -13.476868 22.600006 -1.467890 655.8 163.1

T = length(y) yest = fitted(m) sse = sum((yest-y)ˆ2) ssr = sum((y-mean(y))ˆ2) mse = sse/(T-k) R2 = 1 - sse/ssr Ra2 = 1 - (T-1)*(1-R2)/(T-k) aic = log((T-k)*exp(2*k/T)*mse/T) bic = log(Tˆ(k/T)*(T-k)*mse/T) M = c(Ra2, mse, aic, bic) names(M) = c("R2-ajus","MSE","logAIC","logBIC") return(M)} M.lin = medidas(mod.lin,yi,2) M.cuad = medidas(mod.cuad,yi,3) M.cub = medidas(mod.cub,yi,4) M.exp = medidas(mod.exp,yi,2) M = cbind(M.lin,M.cuad,M.cub,M.exp) (M)

Pr(>|t|) 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.142793 0.000000 0.000000

R2-ajus MSE logAIC logBIC

M.lin 0.892 2.709 1.001 1.018

M.cuad 0.997 0.078 -2.543 -2.517

M.cub 0.997 0.078 -2.544 -2.509

M.exp 0.989 0.264 -1.327 -1.309

Tabla 2.3: Estad´ısticos de Selecci´on de Modelos Del cuadro 2.3 se puede conclu´ır que el modelo que mejor ajusta los datos es el modelo cuadr´atico (2.4) ya que presenta el menor BIC.

Chequeo de las Hip´otesis del Modelo de Regresi´on En este momento el diagn´ostico de incorrelaci´on y normalidad de los errores se hace mediante algunos procedimientos gr´aficos, que se programan con las instrucciones siguientes r = mod.cuad$residuals par(mfrow=c(2,2)) plot(t,r,type=’o’,ylab=’residuo’) abline(h=0,lty=2) plot(density(r),xlab=’x’,main= ’’) qqnorm(r) qqline(r,col=2) acf(r,ci.type="ma",60)

C´alculo de Pron´osticos con el Modelo Cuadr´atico En esta u´ ltima parte se examina c´omo es la calidad de los pron´osticos con el modelo cuadr´atico (2.4). Se utiliza la t´ecnica de validaci´on cruzada que consiste en calcular los pron´osticos con el modelo ajustado con el primer grupo de datos y compararlos con los datos reales del segundo grupo. Se hace con las siguientes instrucciones. # Pron´ osticos tt = seq(482,493,1) tt2 = ttˆ2

30 Year Unemployed Year Unemployed

1950 1869 1988 2242

1960 271 1989 1869

1970 149 1990 1883

1975 1074 1991 1689

1980 889 1992 1808

1985 2304 1993 2270

Tabla 2.4: Unemployed2 Data. pr2 = predict(mod.cuad,data.frame(t=tt,t2=tt2)) plot(tt,yf,type="b") lines(tt,pr2,col="red")

188000 182000

184000

186000

190000

192000

Los pron´osticos est´an en la variable “pr2” generada en la u´ ltima parte del programa. A partir del examen de la figura (2.3) se puede conclu´ır que el modelo cuadr´atico genera pron´osticos confiables a corto plazo.

482

484

486

488

490

492

Figura 2.3: Pron´osticos (—-) versus Observados (-o-o-)

2.5.

Problemas

1. (Unemployed2 Data) Table 2.4 lists total numbers of unemployed (in thousands) in West Germany between 1950 and 1993. Compare a logistic trend function with an allometric one. Which one gives the better fit? 2. Consumo de agua mensual en el a´ rea metropolitana en metros c´ubicos, por estrato, para los estratos 1,2,3,4,5,6, incluyendo el consumo total. Desde enero de 2002. En

31 total son 7 series que se identificar´an con los nombres : consumo 1, consumo 2,..., consumo 6, consumo total. Los datos est´an en el archivo consumo.dat, en la p´agina del curso. Dise˜ne un plan inicial de an´alisis y descr´ıbalo la introducci´on del trabajo. Por ejemplo, analizar las tendencias de todos los estratos, determinar si hay estacionalidad en todos o algunos. a) Lea los datos a partir del archivo asignado que se encuentra en la p´agina web del curso. Grafique la serie. Desarrolle los siguientes puntos y reporte los resultados de cada uno. Con estos resultados elabore el reporte del trabajo. b) Use la funci´on de descomposici´on stl para analizar lo siguiente: 1) Hay una tendencia global que puede describirse por una funci´on lineal, cuadr´atica o´ c´ubica? 2) La tendencia se puede describir mejor por secciones? Es decir, no hay tendencia global? 3) Se puede observar una componente estacional?. Cu´al es el per´ıodo?. c) Si observa que es factible que existan tendencia y estacionalidad, est´ımelas conjuntamente, con la restricci´on que se defini´on para la estacionalidad de considerar solamente s − 1 componentes estacionales (variables indicadoras), donde s es el per´ıodo estacional. Reporte una tabla con los par´ametros estimados, errores est´andar, estad´ısticos t y valores p. Reporte el AIC y el R-cuadrado ajustado. d) Si observa solamente tendencia o´ estacionalidad, est´ımelas y reporte los estad´ısticos del punto anterior. En el caso de estacionalidad reporte el gr´afico el patr´on estacional, es decir, los coeficientes δbi versus i = 1, 2, . . ., s. e) Reporte qqplot normal e histograma de los residuos, as´ı como la prueba JarqueBera de normalidad. Comente sobre estos resultados.

f ) Calcule pron´osticos con base en el modelo estimado. Para esto escoja un horizonte de pron´osticos adecuado a su serie. g) Reporte las conclusiones del trabajo: un resumen de los resultados que encontr´o, y de los problemas que se presentaron en la elaboraci´on del modelo. Por ejemplo, un comentario acerca de lo que Usted crea que logr´o el modelo: captur´o la din´amica de la serie?, su tendencia?, los pron´osticos parecen realistas y confiables?. Qu´e otras alternativas podr´ıan haberse propuesto?. Cu´al es una de las cr´ıticas al modelo de descomposici´on?.

CAP´ITULO

Recursos en R Para Analisis ´ de Tendencia

3.1. Introducci´on En este cap´ıtulo se introducen algunas de las herramientas de R para la implementaci´on de la descomposici´on Yt = Tt + St + εt , y para la estimaci´on de la tendencia Tt. Hay que hacer e´ nfasis en que estos m´etodos no asumen un modelo global para la tendencia Tt sino local, es decir, no son modelos con par´ametros fijos, de la forma, por ejemplo, Tt = β0 + β1 t, sino que β0 y beta1 cambian en el tiempo para permitir mayor flexibilidad. 1. La funci´on loess() realiza un suavizamiento con base en una regresi´on local Loess. 2. Las funciones decompose() y stl() ambas realizan una descomposici´on de la serie Yt en las tres componentes. 3. Hay librer´ıas con varios tipos de filtros para suavizamiento y extracci´on de componentes de tendencia y estacionales como mFilter, que tiene implementado el filtro Hodrick-Prescott. 4. Tambi´en la librer´ıa timsac tiene la funci´on decomp() que realiza una descomposici´on incluyendo una componente autoregresiva y otra para fechas de intervenciones, Yt = Tt + St + Rt + T At + εt . 33

34 5. Mediante la funci´on filter() de la librer´ıa stats se pueden implementar varios tipos de medias m´oviles y filtros lineales, por ejemplo, fitros tipo Henderson y filtros recursivos, que son herramientas u´ tiles para la estimaci´on de la tendencia Tt . En la secci´on 2.1 se introducen los m´etodos de descomposici´on con base en suavizadores . En la secci´on definici´on 2.2 se definen las medias m´oviles. EL uso de las funciones filter (stats) y filter (signal) se describe en la secci´on 2.3.

3.2.

Metodologias de Descomposici´on con Base en Suavizadores

3.2.1.

Regresi´on Local Loess

La Regresi´on Loess (inicialmente Lowess: Locally wighted scatterplot smoothing) es un modelo de regresi´on no-param´etrica que regresa Yi versus xi , pero no asume un modelo global fijo, es decir, no asume un intercepto y una pendiente fijas sino variables, de manera local. Aqu´ı local significa una ventana m´ovil que contiene un n´umero determinado de datos de la variable independiente. La descripci´on del algoritmo Loess es la siguiente. Suponga (xi, Yi ), i = 1, . . ., n. El objetivo de Loess es calcular una funci´on de regresi´on local g(x) que es el an´alogo de, por ejemplo, a + bx, es decir, Yi = g(xi) + εi , i = 1, . . . , n. Suponga x(1) ≤ x ≤ x(n) , entonces g(x) se calcula as´ı: 1. Se escoge q ∈ Z+ tal que q ≤ n. 2. Se escogen los q valores xi m´as cercanos a x  (1 − x3 )3 , 0 ≤ x ≤ 1 3. Defina w(x) = 0 , x≥1

4. Defina λq (x) la distancia de x al xi m´as alejado entre los q escogidos. |xi − x| 5. Defina vi (x) = w , para cada xi escogido. λq (x) 6. Ajuste Yi = a + bxi o´ Yi = a + bxi + cx2i con MCP ponderando cada xi con vi (x). 7. Defina g(x) como el valor estimado de la regresi´on en x. En series de tiempo se toma, por ejemplo, (Yi , ti ), es decir, xi = ti . En cada t se reemplaza Yt por Yˆt = gˆ(t).

100

200

300

400

500

600

700

Ejemplo 3.2.1. Si se toman los datos del n´umero de licencias para vivienda nueva, mensuales, en Medell´ın entre 01/1986 y 06/2003, se observa un fuerte aumento durante los a˜nos 1986 - 1993. Los datos se tomaron de la p´agina web de Camacol (1 ). A partir de 1993 empieza un decrecimiento. Sin embargo, a partir de 1996 se observa un aparente cambio de concavidad y se sugiere una recuperaci´on.

12/85

08/87

04/89

12/90

08/92

04/94

12/95

08/97

04/99

12/00

08/02

fechas

Figura 3.1: Suavizamiento Loess de N´umero de Licencias de Construcci´on en Medell´ın, 1986-2003 Este es el c´odigo en R para suavizamiento con Loess, ver Figura 3.1. G = read.table("licenciasmedellin.dat", header = TRUE, stringsAsFactors = FALSE) attach(G) yw=loess(y ˜ time(y)) np = length(y) fecha = seq(as.Date("1986/01/01"), as.Date("2003/06/03"), by="months") fechas = strptime(as.character(fecha), "%Y-%m-%d") plot(fechas,y, xaxt="n",panel.first = grid(),type=’l’,ylab=’’) axis.POSIXct(1, at=seq(as.Date(fechas[1]),as.Date(fechas[np]), "months"), format="%m/%y") axis.POSIXct(1, at=seq(as.Date(fechas[1]),as.Date(fechas[np]), "years"), labels = FALSE, tcl = -0.2) 1

www.camacol.com.co

36 lines(fechas,yw$fitted, xaxt="n", panel.first = grid(), type=’l’,col=’red’,lwd=2)

3.2.2.

STL: M´etodo de Descomposici´on Basada en Regresi´on Loess

STL es un m´etodo para estimar las componentes Tt y St con base en Regresi´on Loess, desarrollado por Cleveland et al. [1990]. STL consiste de una secuencia de dos aplicaciones iteradas de Regresi´on Loess. Para aplicar este m´etodo se debe especificar una frecuencia de muestreo relacionada con el per´ıodo de la componente estacional. La forma de especificar esta frecuencia es declarando la variable en donde se encuentran los datos como un objeto ts con frecuencia (52, 12, 4, 1), es decir, semanal, mensual, trimestral o anual, respectivamente.

100

200

300

400

500

600

700

Ejemplo 3.2.2. Al aplicar el m´etodo Stl a los datos del Ejemplo (3.2.1) y graficar la serie suavizada se obtiene el resultado de la Figura (3.2).

12/85

08/87

04/89

12/90

08/92

04/94

12/95

08/97

04/99

12/00

08/02

fechas

Figura 3.2: Suavizamiento STL de N´umero de Licencias de Construcci´on en Medell´ın, 1986-2003

3.3.

Filtros Lineales, Medias M´oviles y Suavizadores

Definici´on 3.3.1 (Filtro Lineal). Un filtro lineal es un operador que convierte una serie Yt en otra Y¯t tal que m X Y¯t = wj Yt−j , t = −m, . . . , T − m (3.1) j=−m

37 donde w−m , w−m+1 , . . . , w0, . . . , wm son pesos predeterminados. En algunos casos para determinar la tendencia y reducir las fluctuaciones locales (suavizar) se restringen P a m ovil. Usualmente se escogen r=−m wr = 1. Este filtro lineal se denomina Media M´ w−j = wj , pesos sim´etricos. Los filtros lineales se utilizar´an como estimadores de la tendencia, es decir, Tbt = Y¯t , cuando se observe que la tendencia no parece factible que sea descrita por una u´ nica funci´on P param´etrica, por ejemplo de la forma polin´omica Tt = β0 + kj=1 βj tj .

Una aplicaci´on del filtro lineal es remover la tendencia, es decir, si Tˆt = Y¯t , la serie Xt = Yt − Tˆt es la serie sin tendencia. Si se aplica (3.1) a Yt = Tt + St + εt el filtro lineal altera la amplitud de las oscilaciones de St. Ese efecto se llama ganancia del filtro. Ejemplo 3.3.1. En Diebold [1999, p´ag. 84], se define: La Media M´ovil General Bilateral con wj = w−j , j = 0, 1, . . ., m. Y¯t =

m X

wj Yt−j .

(3.2)

j=−m

El t´ermino Yt−j se llama “rezago de orden j”. Si j < 0 el rezago Yt−j es un valor “futuro”. La Media M´ovil Bilateral wj =

1 , j = −m, . . . , m. 2m + 1

(3.3)

1 , j = 0, 1, . . ., m. m+1

(3.4)

La Media M´ovil Unilateral wj =

El entero m ≥ 1 es el “ancho de ventana” de la media m´ovil. La Media M´ovil Unilateral (3.4) utiliza valores presentes y pasados de la serie. Es un ejemplo de filtro lineal “causal” . La Media M´ovil Bilateral (3.3) no es un filtro lineal causal. Ejemplo 3.3.2. Un ejemplo de filtro lineal es el Suavizamiento Exponencial Simple (SES), definido por m X ¯ Yt = α(1 − α)j Yt−j , α ∈ (0, 1), m > 1. (3.5) j=0

El m´etodo SES tambi´en se denomina con las siglas EWMA (Exponentially Weighted Moving Average). En R se coloca m = b 2−α α c.

38 La ecuaci´on (3.5) se puede obtener de la relaci´on recursiva Y¯t = αYt + (1 − α)Y¯t−1 para P j t = 2, 3, . . . , T . Se comprueba que Y¯t = (1 − α)m+1 Y¯t−m−1 + α m j=0 (1 − α) Yt−j . Si m P j es grande se puede aproximar (1 − α)m+1 ≈ 0 y por tanto Y¯t = m j=0 α(1 − α) Yt−j . Si α ≈ 0, entonces el filtro reduce la serie casi a una constante. De la relaci´on P j ¯ ¯ Y¯t = (1 − α)m+1 Y¯t−m−1 + α m j=0 (1 − α) Yt−j , si α = 0, entonces Yt = Yt−m−1 para todo m > 1. Con lo cual Y¯t = Y¯0 .

Si α ≈ 1, entonces Y¯t = Yt , es decir, se reproduce la serie y no la suaviza. Con α = 0.2 ya se observa esta caracter´ıstica. En R la funci´on emaTA(x, lambda = 0.019, startup = 0) de la librer´ıa fTrading P 2 j , y a la funci´on implementa EWMA Y¯t = m ormula α = m+1 j=0 α(1 − α) Yt−j . Se usa la f´ se le puede ingresar m en lugar de α. El siguiente es un ejemplo de c´odigo en R para EWMA. Ejemplo 3.3.3. C´odigo en R para SES (EWMA) con emaTA. Los resultados est´an en la Figura 3.3. N´otese que en “lambda” se ingresa el valor de α.

library(fTrading) x = MSFT # Base de Datos de Microsoft OHLC x = x[,"Close"] # Precios de Cierre y = emaTA(x, lambda = 0.189) seriesPlot(x) lines(y, col="red")

2000−09−27

2001−02−20

(a) α = 0.019

2001−07−16

2000−09−27

2001−02−20

2001−07−16

(b) α = 0.189

Figura 3.3: Precios de la Acci´on de Microsoft Suavizados con SES

3.4. La Funci´on Filter Una funci´on u´ til para programa medias m´oviles y suavizadores es la funci´on filter. En esta secci´on se explica su funcionamiento y aplicaciones. En R existen dos versiones de filter: La primera est´a en la librer´ıa stats que es cargada por defecto en R, y la segunda est´a en la librer´ıa signal. La exposici´on siguiente es sobre la funci´on filter de la librer´ıa stats. Cuando se utilice la de la librer´ıa signal se escribe signal::filter(). Esta u´ ltima se describe en la secci´on de Notas. La forma de utilizar esta funci´on es la siguiente. x = filter(y, w, method = c("convolution", "recursive"), sides = 1,2, circular = TRUE, FALSE, init) Argumentos

y: w: method:

sides: circular: init:

Es la serie que se quiere suavizar.Puede ser multivariada. Es el vector con los pesos w = (w0 , w1, . . . , wm), con wj real o´ matriz. Cualquiera de las dos "convolution" o "recursive". Si es "convolution" la funci´on calcula una media m´ovil. Si es "recursive" calcula un filtro recursivo. P Si se usa "convolution", y sides=1 calcula m j=0 wj Yt−j Pm Si sides=2 calcula j=−m w|j| Yt−j . Solo cuando se usa "convolution", es TRUE o´ FALSE. Solo cuando se usa "recursive".

La opci´on method = “convolution” La opci´on method = "convolution" se utiliza para calcular medias m´oviles. 1 Pm Ejemplo 3.4.1. Para calcular Y¯t = m+1 j=0 Yt−j se programa: x = filter(y,rep(1/(m+1),m+1),"conv",1,T,NULL) P Ejemplo 3.4.2. Para calcular Y¯t = m j=0 wj Yt−j

w = c(w0, w1,...,wm) x = filter(y, w, "conv",1,T,NULL)

40 La opci´on "circular=TRUE" requiere m´as explicaci´on. En la media m´ovil unilateral P se calcula Y¯t = m j=0 wj Yt−j . Si se escribe la serie en la forma Yt−T +1 , . . . , Yt−1 , Yt , que es como la lee filter, y se tienen los m + 1 pesos w0 , w1 , . . . , wm, al aplicar X=filter(Y,w,"conv",1,T,NULL) se obtiene una serie X = (Xt−T +1 , . . . , Xt−1 , Xt), calculada de la siguiente manera. Xt = w0 Yt + w1 Yt−1 + . . . + wmYt−m Xt−1 = w0 Yt−1 + w2 Yt−2 + w3 Yt−3 + . . . + wm Yt−m−1 .. . Xt−T +1 = w0 Yt−T +1 + w1 Yt−T + w2 Yt−T −1 + . . . + wm Yt−T −m+1

Pero los datos Yt−T , . . . , Yt−T −m+1 no est´an disponibles. El efecto de la opci´on TRUE para "circular" es reemplazarlos por los valores Yt , Yt−1 , . . ., Yt−m+1 , respectivamente. De manera similar se reemplazan en las medias bilaterales. Ejemplo 3.4.3. Ejemplo simple de convolution con reemplazo de valores. y = c(1,1,3,1); w = c(1,2,3) x = filter(y,w,"conv",1,T,NULL) # x tiene cuatro elementos calculados as´ ı x[4] = 1(1) + 2(3) + 3(1) = 10 x[3] = 1(3) + 2(1) + 3(1) = 8 x[2] = 1(1) + 2(1) + 3(1) = 6 x[1] = 1(1) + 2(1) + 3(3) = 12 P Ejemplo 3.4.4. Para calcular Y¯t = m j=−m w|j| Yt−j

# k = 2m+1 w = c(w0,w1, w2,...,wk) x = filter(y,w,"conv",2,T,NULL)

Ejemplo 3.4.5. C´alculo de EWMA con la funci´on filter. Como en EWMA (3.5) se tiene Xt = αYt + α(1 − α)Yt−1 + α(1 − α)2 Yt−2 + · · · + α(1 − α)m Yt−m , se puede generar X as´ı: m = floor(1/a - 1) w = a*(1-a)ˆseq(0,m,1) x = filter(y,w,"conv",1,T,NULL)

41 Ejemplo 3.4.6. La siguiente media m´ovil de orden 12 estima la componente Tt en una serie con componente estacional de per´ıodo p = 12. 5 X 1 1 1 Yt−6 + Tt = Yt−u + Yt+6 , 12 2 2

t = 7, . . . , 45,

(3.6)

u=−5

Aplicar este filtro a la serie de tiempo “nottem”, que est´a por defecto en R, y corresponde a las temperaturas promedio mensuales, en grados Farenheit, del castillo de Nottingham, en Inglaterra. La programaci´on en R del filtro (3.6) se puede hacer de la manera siguiente. #--------------------------Yt = nottem ts.plot(Yt) f = frequency(Yt) w = c(0.5, rep(1, f - 1), 0.5)/f Tt = filter(Yt,w,"conv",2,T,NULL) Xt = Yt - Tt t=seq(1,length(Yt)) par(mfrow = c(3,1)) plot(t,Yt,type=’l’); plot(t,Tt,type=’l’); plot(t,Xt,type=’l’); # Notese que el filtro estima la componente de tendencia # exactamente como lo hace la funci´ on stl() m = stl(Yt, "per") Tmt = m$time.series[,2] par(mfrow = c(1,1)) plot(t,Tt,type=’l’); lines(t,Tmt,col=’red’);

La opci´on method= “recursive” La opci´on method = "recursive" calcula una serie X definida por Xt =

m X j=1

wj Xt−j + Yt ,

(3.7)

42 donde los pesos wj se deben ingresar en un vector de la forma (w1, . . . , wm).A partir de la serie Y = (Yt−T +1 , YtT , . . ., Yt−1 , Yt) calcula X = (Xt−T +1 , . . . , Xt−1, Xt). N´otese que P Xt−T +1 = m j=1 wj Xt−j−T +1 + Yt−T +1 requiere los valores Xt−j−T +1 , j = 1, . . ., m. Pero estos m valores no est´an disponibles y deben proporcionarse en un vector de valores iniciales, por ejemplo Z = (Z1 , . . . , Zm ). El programa se muestra a continuaci´on.

Ejemplo 3.4.7. Para calcular Xt =

j=1 wj Xt−j

+ Yt

w = c(w1,...,wm) z = c(z1,...,zm) x = filter(y, w, "rec",init=z)

Ejemplo 3.4.8. Para calcular Xt = de “conv” y “rec” como sigue.

a b z u x

= = = = =

3.5.

j=1

aj Xt−j +

j=0 bj Yt−j

se hace una combinaci´on

c(a1,...,ap) c(b0,...,bq) c(z1,...,zp) filter(y, b, "conv",1,T,NULL) filter(u, a, "rec", init=z)

Ejemplo de An´alisis de Tendencia con Suavizadores

Este programa calcula varios suavizamientos de la serie del ´ındice de la Bolsa de Colombia, IGBC, durante el per´ıodo 19/01/2009 a 19/08/2010, que fu´e un per´ıodo al alza. Los resultados de muestran a en la Figura 3.4 continuaci´on

obs ses holt m.doble

01/09

03/09

05/09

07/09

09/09

11/09

01/10

03/10

05/10

07/10

fechas

Figura 3.4: Suavizamientos SES, Holt, Media Movil unilateral 30 del IGBC, 2009-2010 # analisis del igbc mediante ewma, medias moviles # bi-unilaterales y filter library(forecast) library(fTrading) source("holt.r") D = read.table("igbc.dat", header = TRUE, stringsAsFactors = FALSE) attach(D) n1 = which(fecha=="19/01/2009") n2 = which(fecha=="19/08/2010") x = igbc[n1:n2]/1000 f = fecha[n1:n2] # suavizamiento exponencial EWMA f22 = emaTA(x, lambda = 0.072, startup = 30) # suavizamiento Holt f.obj = function(theta,yt){

44 a = theta[1] b = theta[2] n = length(yt) yet = holt(yt,a,b,n)$yet err = yt[2:n]-yet[1:(n-1)] mse = mean(errˆ2) return(mse) } # parametros que minimizan el MSE (alpha = optim( c(0.5,0.5), f.obj, lower = rep(1.0e-06,2), upper = rep(0.9999999,2), method = "L-BFGS-B", yt=x)$par) # otro par de parametros alfa1 = 0.032 alfa2 = 0.0003 np = length(x) f33 = holt(x,alfa1,alfa2,np)$yet # medias bilaterales (filtro henderson) w13 = c(-0.019, -0.028, 0.0, 0.066, 0.147, 0.214, 0.240, 0.214, 0.147, 0.066, 0.0, -0.028, -0.019) w23 = c(-0.004, -0.011, -0.016, -0.015, -0.005, 0.013, 0.039, 0.068, 0.097, 0.122, 0.138, 0.148, 0.138, 0.122, 0.097, 0.068, 0.039, 0.013, -0.005, -0.015, -0.016, -0.011, -0.004) fw13 = filter(x,w13,"conv",2,F,NULL) fw23 = filter(x,w23,"conv",2,F,NULL) f4 = filter(x,rep(1/30,30),"conv",2,F,NULL)

np = length(x) fecha = as.Date(f, "%d/%m/%Y") fechas = strptime(as.character(fecha), "%Y-%m-%d") plot(fechas,x, xaxt="n",panel.first = grid(), type=’l’,ylab=’’,col=’grey’) axis.POSIXct(1, at=seq(as.Date(fechas[1]), as.Date(fechas[np]), "months"), format="%m/%y") axis.POSIXct(1, at=seq(as.Date(fechas[1]), as.Date(fechas[np]), "years"), labels = FALSE, tcl = -0.2) lines(fechas,f22, xaxt="n", panel.first = grid(), type=’l’,lty=3,lwd=2,col=’darkgrey’) lines(fechas,f33, xaxt="n", panel.first = grid(), type=’l’,lty=5,lwd=2,col=’darkgrey’) lines(fechas,f4, xaxt="n", panel.first = grid(), type=’l’,lty=7,lwd=2,col=’darkgrey’) legend("topleft", c("obs","ses","holt", "m.doble"), lty = c(1,3,5,7), lwd=c(1,2,2,2) )

plot(fechas,x, xaxt="n",panel.first = grid(), type=’l’,ylab=’’,col=’grey’) axis.POSIXct(1, at=seq(as.Date(fechas[1]), as.Date(fechas[np]), "months"), format="%m/%y") axis.POSIXct(1, at=seq(as.Date(fechas[1]), as.Date(fechas[np]), "years"), labels = FALSE, tcl = -0.2)

lines(fechas,f13, xaxt="n", panel.first = grid(), type=’l’,lty=3,lwd=2,col=’darkgrey’) lines(fechas,f23, xaxt="n", panel.first = grid(), type=’l’,lty=5,lwd=2,col=’darkgrey’) lines(fechas,f4, xaxt="n", panel.first = grid(),

46 type=’l’,lty=7,lwd=2,col=’darkgrey’) legend("topleft", c("obs","hender13","hender23", "m.doble"), lty = c(1,3,5,7), lwd=c(1,2,2,2) )

3.6.

Notas

Sobre la funci´on filter de la librer´ıa signal La librer´ıa signal tiene otra versi´on de la funci´on filter. La forma de aplicar filter es x = signal::filter(b,a,y). Los vectores a y b se definen como: a = c(1, a1, . . . , ap) b = c(b0 , b1, . . . , bq ) El vector a se denomina de coeficientes recursivos, el vector b se denomina de coeficientes de media m´ovil. Note que la primera componente del vector a es 1. Al aplicar x=signal::filter(b,a,y) a Y = (Yt−T +1 , YtT , . . ., Yt−1 , Yt) se obtiene una serie filtrada X = (Xt−T +1 , . . . , Xt−1 , Xt) definida por: Xt +

p X

aj Xt−j =

j=1

q X

bj Yt−j .

(3.8)

j=0

P P Ejemplo 3.6.1. Suponga que se quiere calcular Xt = pj=1 aj Xt−j + qj=0 bj Yt−j . Se puede utilizar la funci´on filter de la librer´ıa signal o´ la funci´on filter de la librer´ıa stat. Con la librer´ıa stat el programa es el mismo del ejemplo (3.4.8). La programaci´on con la funci´on filter de signal es como sigue. a b z x

= = = =

c(a1,...,ap) c(b0,...,bq) c(z1,...,zp) signal::filter(b, c(1,-a),y)

N´otese que se cambia el signo al vector a para que al reemplazar en la definici´on en (3.8) se obtenga el resultado deseado. Adem´as, n´otese que al cargar la librer´ıa signal se debe diferenciar la funci´on filter indicando de cu´al librer´ıa se trata colocando x=signal::filter. No parece que haya una ventaja muy grande al utilizar esta versi´on de filter.

47 Sobre EWMA N´otese que tambi´en a partir de Xt = αYt + (1 − α)Xt−1 , se puede programar w x n m x

= = = = =

1-a filter(a*y,w,"rec") length(x) floor((2-a)/a)) x[(m+1):n]

M´etodo de Suavizamiento de Holt El m´etodo de suavizamiento de Holt consiste en un par de ecuaciones recursivas, con un par de par´ametros, (α, β) ∈ (0, 1) × (0, 1). at = αYt + (1 − α)(at−1 + bt−1 ), bt = β(at − at−1 ) + (1 − β)bt−1 .

(3.9) (3.10)

con Ybt = at + bt , para t = 1, . . . , T . Los valores iniciales a0 , b0 se toman como el intercepto y la pendientes estimadas por regresi´on en el modelo Yt = a + bt + εt . Un programa en R para una funci´on en R que retorna (Ybt, at, bt), se muestra a continuaci´on.

#suavizamiento de Holt -------------------holt = function(yt,a,b,n){ bst = double(n) ast = double(n) yet = double(n) rw = lm(yt ˜ time(yt)) ast[1] = rw$coef[1] bst[1] = rw$coef[2] yet[1] for (i { ast[i] bst[i] yet[i]

=yt[1] in 2:n) = a*yt[i]+(1-a)*(ast[i-1]+bst[i-1]) = b*(ast[i]-ast[i-1])+(1-b)*bst[i-1] = ast[i]+bst[i]

48 } D = data.frame(yet,ast,bst) return(D)}

CAP´ITULO

Modelado y Pronostico ´ de Series Estacionales

4.1. Definiciones En este cap´ıtulo se introducen modelos para la componente estacional St que forma la descomposici´on de una serie Yt = Tt + St + εt , con Tt la tendencia y εt la componente aleatoria. Adem´as se exponen procedimientos de estimaci´on de St con base en regresi´on lineal y filtros. Definici´on 4.1.1 (Componente Estacional). La componente estacional St se define como una funci´on no aleatoria, peri´odica de per´ıodo s. Los valores de St para t = 1, . . ., s se denominan “el patr´on estacional”. El per´ıodo estacional s es el n´umero m´ınimo de per´ıodos que tarda el patr´on estacional en repertirse. La unidad de tiempo t de Yt es d´ıa, semana, mes, trimestre, semestre, a˜no. Definir el per´ıodo estacional por ejemplo como s = 12 significa que la unidad de tiempo es mes, y se muestrea doce veces al a˜no. De manera similar con s = 4, trimestral. Un per´ıodo semanal es s = 52 y se muestrea 52 veces al a˜no. Tambi´en se puede definir un modelo con t en d´ıas, y s = 7, en el cual se muestrea 360 = 52(7), veces al a˜no y se define un patr´on estacional que se repite cada 7 d´ıas. Ejemplo 4.1.1. Una ejemplo de series con componente estacional es la serie nottem, en la librer´ıa stats de R. Corresponden a la serie mensual “Average Monthly Temperatures 49

50 30

nottem

at Nottingham, 1920-1939”, que contiene las temperaturas promedio de grados Fahrenheit tomadas en el Castillo de Nottingham (UK), por 20 a˜nos, entre 1920 y 1939. Datos del libro de Anderson, O. D. (1976) Time Series Analysis and Forecasting: The Box-Jenkins Approach. La componente estacional de per´ıodo s = 12 se observa claramente. (1 )

1920

1925

1930

1935

1940

Time

Figura 4.1: Average Monthly Temperatures at Nottingham, 1920-1939

Propiedades de St En todo lo que sigue St , t = 1, 2, . . . es una funci´on real en tiempo discreto, es decir, es una sucesi´on. Que cumple las propiedades siguientes. 1. St es una funci´on peri´odica con periodo s, St+s = St para t = 1, 2, . . .. Por tanto, s´olo se requiere definir St en los primeros s valores, St , t = 1, 2, . . . , s. Es decir, basta con definir el patr´on estacional. 2. Si S1t y S2t son funciones estacionales con periodo s entonces aS1t + bS2t, para a, b ∈ R, es tambi´en una funci´on estacional de per´ıodo s. El conjunto de funciones peri´odicas de per´ıodo s es un espacio vectorial de dimensi´on s. 3. Una base para este espacio vectorial est´a conformada por las s variables indicadoras estacionales, Ij (t) ∈ {0, 1}, para j = 1, . . . , s, t = 1, 2, . . ., definidas por ( 1 , t = j, j + s, j + 2s, . . . (4.1) Ij (t) = 0 , en otro caso 1

http://www.astrostatistics.psu.edu/datasets/R/html/datasets/html/00Index.html

51 O, de manera equivalente, Ij (t) = 1, si t = sbt/sc + j, donde bxc es la parte entera de x, Ij (t) = 0 en otro caso. Por ejemplo, si s = 4, j = 3, t = 11, entonces como 11 = 4b11/4c + 3, se tiene I3 (11) = 1. Tambi´en se puede usar el concepto de “congruencia m´odulo s”, Ij (t) = 1 si t ≡ j(mod s), es decir, si t − j es divisible por s. Por tanto, si St es peri´odica de per´ıodo s, existen constantes δj , j = 1, . . . , s tales que s X St = δj Ij (t), t = 1, 2, . . .. (4.2) j=1

4. Tambi´en se puede utilizar el conjunto {sen(2πjt/s), cos(2πjt/s), j = 1, . . ., s}. Para St peri´odica de per´ıodo s, escogemos m ≥ 1 y constantes b1,j , b2,j , j = 1, . . ., m tales P que, St ≈ m j=1 β1,j sin(2πjt/s) + β2,j cos(2πjt/s), para t ≥ 1. El valor de m se escoge 1 ≤ m ≤ 2s − 1, si s es par, o por 1 ≤ m ≤ s−1 2 si s es impar. Valores m = 1, 2 son recomendables por parcimonia. Note que con m = 2 habr´ıan cuatro par´ametros β1,j , β2,j a estimar mediante regresi´on lineal. 5. N´otese que en el modelo con variables indicadoras el n´umero de par´ametros aumenta con el per´ıodo s. Para per´ıodos grandes, por ejemplo s = 360 no ser´ıa un modelo pr´actico y ser´ıa preferible el modelo con variables trigonom´etricas. P Si se define un modelo para la componente de tendencia Tt, por ejemplo, Tt = kj=0 βj tj , k = 1, 2, 3, y el modelo para la componente estacional con variables indicadoras (4.2), entonces el modelo para Yt est´a dado por: Yt = β0 +

k X j=1

βj t +

s X

δj Ij (t) + εt ,

(4.3)

j=1

donde k = 1, 2, 3 y βj , j = 1, . . ., k y δj , j = 1, . . ., s son par´ametros a estimar. Es decir, P se puede escribir Yt de la forma Yt = β0 + pi=1 βi Xi,t + εt , donde Xi,t es de la forma ti o´ Ii (t), y se pueden utilizar todos los recursos de regresi´on lineal m´ultiple para la estimaci´on y pron´ostico con este modelo.

4.2. Estimaci´on de la Componente Estacional 4.2.1.

Modelos con Variables Indicadoras

A continuaci´on se enumeran los pasos para estimar la componente St, asumiendo primero el modelo con variables indicadoras (4.2).

52 1. Identificar el per´ıodo de la serie Yt , s, si es que existe. Con R al declarar la serie como objeto “ts” con frecuencia, por ejemplo 12, ya est´a idenficado el posible per´ıodo como s = 12. 2. Generar las s variables indicadoras Ij (t) para t = 1, 2, . . ., T . Es decir, generar la matriz Ij (t), donde al variar j = 1, . . . , s variamos las columnas y al variar t = 1, . . ., T varian las filas. P 3. Estimar los par´ametros en el modelo Yt = Tt + sj=1 δj Ij (t) + εt mediante regresi´on lineal. N´otese que se trata de estimar los par´ametros de las componentes de tendencia y estacional, conjuntamente. La funci´on Tt + St se denominar´a “componente estructural” , y los residuales εt se denominar´an “residuos estructurales” .

El Problema de la Multicolinealidad con Variables Indicadoras Se deben inclu´ır solamente s − 1 variable indicadoras en el modelo en lugar del conjunto completo de s, para evitar problemas de multicolinealidad ya que las variables indicadoras P cumplen sj=1 Ij (t) = 1, ∀t, por lo que estar´ıan perfectamente correlacionadas con el intercepto β0 , pues, aunque no se hace expl´ıcito, esta constante est´a asociada con una variable explicativa que es constante e igual a 1. El modelo(4.3) se modifica y queda Yt = β0 +

k X j=1

βj tj +

s−1 X

δj Ij (t) + εt ,

(4.4)

j=1

ver Diebold [1999, p´ag. 90].

Estimaci´on en R del modelo con Variables Indicadoras En la librer´ıa forecast hay varios recursos y funciones para estimar y pronosticar con el modelo (4.3). Suponga que Yt tiene una componente de per´ıodo s que se declara con frecuencia s. El primer paso consiste en generar las variables indicadoras Ij(t), j = 1, . . . , s, t = 1, . . ., T . La funci´on seasonaldummy de la librer´ıa forecast calcula una matriz de ceros y unos, denotada por ejemplo, It, de dimensi´on T × (s − 1), con s − 1 columnas, a partir de la

53 informaci´on sobre la serie Yt de frecuencia s. Si el dato de Y1 corresponde a la estaci´on j entonces la primer fila de la matriz It ser´a un vector de ceros excepto en la posici´on j, que corresponde a un uno. De acuerdo a este valor inicial se genera la matriz It. Adicionalmente, el comando seasonaldummyf genera otra matriz de ceros y unos, denotada por ejemplo Itp, con s − 1 columnas y con un n´umero de filas igual al n´umero de pron´osticos a partir de T . El patr´on de las estaciones en esta matriz est´a de acuerdo a la informaci´on de la serie, como se indic´o. Antes de un ejemplo con seasonaldummy es interesante analizar el siguiente c´odigo en R que calcula la matriz Ij (t) de dimensi´on T × s. (2 ). Ejemplo 4.2.1. #-------------------------dummy=function(s,n,i){ # s: estaciones, n: numero de datos, i: estaci´ on de inicio A=diag(s) for(j in 1:(floor(n/s))){ A=rbind(A,diag(s))} A=A[i:(n+i-1),] return(A) } Ejemplo 4.2.2. Ejemplo con R para modelaci´on de la serie de producci´on de cemento Portland, trimestral, en toneladas por mil, entre Q1 1956 y Q3 1994, en Australia. La serie se asume que es de la forma Yt = β0 + β1 t +

s−1 X

δj It (t) + εt .

(4.5)

j=1

Es decir, el modelo a estimar es lineal para la tendencia, m´as una componente estacional. # ----------------------------------------library(forecast) library(xtable) E = read.table("cementq.dat", header = TRUE) attach(E) y = ts(y,frequency=4,start=c(1956,1),end=c(1994,3)) #----modelo con tendencia lineal y estacionalidad #----con variables indicadoras estacionales 2

autor : Ledwing Osorio C´ardenas

t = seq(1,length(y)) It = seasonaldummy(y) mod1 = lm(y ˜ t + It) summary(mod1) print(xtable(mod1),digits=6) r1 = mod1$residuals yhat1 = mod1$fitted.values #-----------------------------

500

1000

1500

2000

P Los resultados del modelo lineal con variables estacionales Yt = β0 +β1 t+ 3j=1 δj Ij (t)+ εt , est´a en el Cuadro 4.1 siguiente. Una comparaci´on entre la serie ajustada y la estimada 2 se puede ver en la Figura (4.2). Con un R-cuadrado ajustado R = 0.8966, el ajuste resulta aceptable. Sin embargo, la componente estacional St parece no ser constante.

100

150

Figura 4.2: Ajuste con Variables Indicadoras

4.2.2.

Modelos con Funciones Trigonom´eticas

En este caso el modelo a estimar es:

Yt = β0 +

m X j=1

βj tj +

k X j=1

β1,j sen(2πjt/s) + β2,j cos(2πjt/s) + εt ,

(4.6)

Tabla 4.1: Resultados del Ajuste del Modelo con Variables Estacionales

(Intercept) t ItQ1 ItQ2 ItQ3

Estimate 665.2149 7.2364 -163.0881 -33.5553 12.0287

Std. Error 24.1428 0.2023 25.6796 25.6788 25.6796

t value 27.55 35.78 -6.35 -1.31 0.47

Pr(>|t|) 0.0000 0.0000 0.0000 0.1933 0.6402

donde el l´ımite k debe estar acotado por: 1 ≤ k ≤ 2s − 1, si s es par, o por 1 ≤ k ≤ s−1 2 si s es impar. El procedimiento de estimaci´on empieza con la generaci´on de las variables explicativas sen(2πit/s), cos(2πit/s). Para calcular estas variables se utiliza la funci´on It.trig = fourier(y,j), donde j = 1, 2, . . ., la cual calcula los harm´onicos sen(2πjt/s), cos(2πjt/s), para cada t = 1, 2, . . ., T . Ejemplo 4.2.3. Un modelo para la demanda diaria de energ´ıa del mercado Nord Pool (pa´ıses escandinavos) es el siguiente. La serie de precios es diaria para un per´ıodo de 30/12/1996-26/04/2000. Se asume que la serie de precios es la suma de cuatro componente: componente lineal (a + bt), componente estacional semanal (st ), componente estacional anual (St ), y componente aleatoria (et ):

Yt = a + bt + st + St + et , st =

6 X

(4.7)

δj Ij (t),

j=1

St = α cos

2πt 365

+ β sen

2πt , 365

En este modelo la componente con per´ıodo anual se modela con funciones trigonom´etricas y la componente con per´ıodo semanal se modela con variables indicadoras. Cuando se aplica la descomposici´on con la funci´on stl() e´ sta estima la componente de tendencia como a + bt + St. El siguiente c´odigo en R implementa el modelo (4.7). # C´ odigo en R -------------------library(forecast)

56 D = read.table("serie.dat",header=T) attach(D) y = ts(y, frequency = 7) n = length(y) t = 1:n St = seasonaldummy(y) x1 = cos(2*pi*t/365) x2 = sin(2*pi*t/365) modelo1 = ln(y ˜ t + St + x1 + x2) summary(modelo1) Ejemplo 4.2.4. Continuando con el Ejemplo (4.2.2), ahora se estima el modelo (5.13), pero utilizando funciones trigonom´etricas. Yt = β0 + β1 t +

2 X

β1,j sin(2πjt/4) + β2,j cos(2πjt/4) + εt .

(4.8)

j=1

Se incluye el c´alculo de pron´osticos a 8 trimestres con los modelos de variables indicadoras y de funciones trigonom´etricas. Los resultados est´an en el Cuadro 4.2. #----modelo con tendencia lineal y estacionalidad #----con variables trigonometricas It.trig = fourier(y,2) mod2 = lm(y ˜ t + It.trig) summary(mod2) print(xtable(mod2),digits=6) r2 = mod2$residuals yhat2 = mod2$fitted.values # pronosticos 8 trimestres # con el modelo con variables indicadoras T = length(yi) Itp = seasonaldummyf(y,8) tp = seq(T+1,T+8,1) prons.vi = predict(mod1, data.frame(t = tp, It=I(Itp)))

# con el modelo con funciones trigonometricas Itp.f = fourierf(y,2,8) prons.vt = predict(mod2, data.frame(t = tp, It.trig=I(Itp.f))) print(xtable(rbind(prons.vi,prons.vt))) #-------------------------------------

Tabla 4.2: Resultados de ocho pron´osticos con Holt-Winters, Indicadoras y Trigonom´etricas Trimestre 1992 Q4 1993 Q1 1993 Q2 1993 Q3 1993 Q4 1994 Q1 1994 Q2 1994 Q3

Obs 1569.00 1450.00 1569.00 1648.00 1777.00 1468.00 1732.00 1962.00

H-W 1572.62 1320.77 1488.50 1570.55 1602.48 1350.63 1518.35 1600.40

Indicadora 1752.06 1599.55 1733.36 1780.69 1781.63 1629.12 1762.93 1810.26

Trigonom 1750.25 1598.36 1738.47 1771.06 1793.91 1632.41 1763.56 1805.11

La conclusi´on que puede obtenerse es que en este caso es indiferente utilizar variables indicadoras o´ funciones trigonom´etricas, dado que los pron´osticos son muy similares. A esta conclusi´on se puede llegar tambi´en observando que los valores de la ra´ız del error cuadr´atico medio, 111.1831, 111.1278, para el modelo de variables indicadoras y para el de funciones trigonom´etricas, respectivamente, son pr´acticamente id´enticos.

4.3. Procedimientos Recursivos de Estimaci´on para Diagnosticar Estabilidad Estructural en Modelos de Pron´osticos ((Con frecuencia, las relaciones comerciales y econ´omicas vari´an con el tiempo. los procedimientos recursivos de estimaci´on nos permiten evaluar y rastrear los par´ametros variables en el tiempo, y, en consecuencia, son u´ tiles para elaborar y evaluar diversos modelos de pron´osticos.)) Diebold [1999, p´ag. 92]

58 Suponga el modelo de regresi´on lineal Yt =

k X

βi Xi,t + εt ,

t = 1, . . ., T,

(4.9)

i=1

i.i.d.

donde εt ∼ N (0, σ 2). Si lo par´ametros βi no fueran constantes sino que cambiarian con el tiempo t, de manera P aleatoria o no, se tendr´ıa Yt = ki=1 βi,tXi,t + εt , un modelo con coeficientes variables.

Note que (4.9) tambi´en se puede escribir como Yt = X 0t β+εt con X t = (X1,t, X2,t, . . . , Xk,t) vector fila en Rk y β ∈ Rk vector de par´ametros. Pero en caso de que los par´ametros var´ıen en el tiempo el modelo quedar´ıa Yt = X 0t β t + εt . P Definici´on 4.3.1 (Estabilidad Estructural). Si Yt = ki=1 βi Xi,t + εt , t = 1, . . ., T , la estabilidad estructural del modelo se define como el supuesto de que los coeficientes βi se mantengan constantes en cada tiempo t. Lo contrario es que alg´un βi cambie con t. La primera afirmaci´on es la hip´otesis nula y la segunda es la alterna. Ejemplo 4.3.1. Algunos ejemplos de cambio estructural en la literatura, debidos a intervenciones en fechas espec´ıficas, son EL nivel del rio Nilo, serie anual, con un cambio en 1898 con la primera represa de Ashwan. El n´umero de muertes en accidentes automovil´ısticos mensuales en Gran Breta˜na. La medida del cintur´on de seguridad en Enero 31 de 1983 cambi´o la serie mediante una intervenci´on que disminuy´o el nivel. En Harvey and Durbin [1986] se hace un an´alisis completo utilizando modelos estructurales. En Box and Tiao [1975, pag. 70] se exponen dos intervenciones en la ciudad de Los Angeles: la apertura del la autopista Golden State Freeway y la entrada en vigencia de la ley 63 que reduc´ıa la proporci´on de hidrocarburos en la gasolina. Ambas intervenciones se esperaba que produjeran un efecto en los niveles de contaminaci´on del aire en los inicios de la d´ecada de 1960, medidos en la serie de ppm horaria de O3 . Dos t´ecnicas para chequear el supuesto de estabilidad estructural son la estimaci´on recursiva de los par´ametros del modelo y las pruebas tipo CUSUM. El objetivo al aplicar estas t´ecnicas es asegurar que el modelo para pron´osticos es estable. En caso de no tener estabildad estructural los pron´osticos con un modelo de regresi´on global podr´ıan no ser confiables. En este caso ser´ıa recomendable utilizar un procedimiento adaptativo como Holt, Holt-Winters, Loess, o´ un modelo de variables estructurales.

59 Estimaci´on Recursiva de los Par´ametros Definici´on 4.3.2 (Estimadores Recursivos). En lugar de utilizar todos los datos para la estimaci´on (t = 1, . . ., T ), si el modelo tiene k par´ametros se utilizan los k primeros datos y se estiman los par´ametros βˆj,k . Luego se utilizan los primeros k + 1, y se estiman los βˆj,k+1 , y as´ı sucesivamente hasta utilizar los T datos. Los par´ametros obtenidos βˆj,t t = k, . . . , T se denominan Estimadores Recursivos.

Para implementar en R la estimaci´on recursiva se utilizan las siguientes instrucciones. Suponga que el n´umero de par´ametros es k y el total de datos es T y que se trata de un modelo Ps−1 con tendencia lineal y estacional, dado por Yt = β0 + β1 t + j=1 δj Ij (t) + t , t = 1, . . . , T . El total de par´ametros es k = 2 + s − 1. Se toman un n´umero inicial de datos k un poco mayor de 2 + s − 1, por ejemplo, k = 2 + s − 1 + 10, para tratar de evitar fluctuaciones en los estimadores por poca informaci´on, y se corre el siguiente programa.

#--------programaci´ on de los estimadores recursivos k = 2+s-1+10 n = T-k parm = mat.or.vec(n,(2+s-1)) for(j in 1:n){ yj = y[1:(k+j)] tj = t[1:(k+j)] Itj = It[1:(k+j),] mod.j = lm(yj ˜ tj + Itj) parm[j,] = t(mod.j$coefficient) } colnames(parm)=c("beta.0","beta.1", "delta.1","delta.2","delta.3") plot.ts(parm) Ejemplo 4.3.2. Con base en el Ejemplo (4.2.2), de la serie de producci´on trimestral de cemento, la gr´afica de los estimadores recursivos de los 5 par´ametros para el modelo con tendencia lineal y variables indicadoras estacionales (5.13), pag. 86, est´an en la Figura 4.5. Lo que puede conclu´ırse de esta figura es que se aprecia un cambio en las trayectorias de los estimadores, por ejemplo en β0 , en β1 y en δ3 . Por tanto, se puede conclu´ır de este diagn´ostico que existe evidencia de inestabilidad estructural en el modelo.

delta.2 delta.3

5 10 15 20 25 30 35 40 −30 −20 −10

600 9 8

beta.1

11500

550

beta.0

650

delta.1

−160 −140 −120 −100 −80 −60

100

120

140

Time

100

120

140

Time

Figura 4.3: Estimadores Recursivos para el Modelo de Producci´on de Cemento Pruebas CUSUM Definici´on 4.3.3 (Residuales Recursivos.). Para cada t = k, k + 1, . . . , T − 1 se calcula P el pron´ostico a un paso, Yˆt = kj=1 βˆj,t Xj,t donde βˆj,t son los estimadores recursivos en t. y se forma el Residual Recursivo εˆt+1,t = Yt+1 − Yˆt ,

(4.10)

Si σ 2 es la varianza de error, se puede probar que

con rt = 1 + X 0t+1

εˆt+1,t ∼ N (0, σ 2rt ),

j=1

X j X 0j

−1

(4.11)

X t+1 .

Los residuales recursivos estandarizados se definen como εˆt+1,t Wt+1,t = √ , σ ˆ 2 = M SE, σ ˆ rt

(4.12)

i.i.d.

y se tiene que Wt+1,t ∼ N (0, 1). Definici´on 4.3.4 (CUSUM). La suma acumulada de los residuos recursivos o´ CUSUM, se define como t X Wi+1,i , t = k, k + 1, . . . , T − 1, (4.13) CU SU Mt = σW j=1

donde

σW =

PT −1 i=k

¯ Wi+1,i − W T −k

¯ = W

T −1 1 X Wi+1,i . T −k i=k

61 Uso: El estad´ıstico CU SU Mt se utiliza para chequear la estabilidad estructural del modelo.

−3

−6

−2

−4

−1

−2

La regi´on de aceptaci´on de la prueba es una banda formada por 2 l´ıneas rectas.

0.0

0.2

0.4

0.6

0.8

1.0

(a) No se rechaza H0

1880

1900

1920

1940

1960

(b) Se rechaza H0

Figura 4.4: Regi´on de Rechazo de la Prueba CU SU M . La expresi´on para las rectas es ±a

√

t−k T − k + 2√ T −k

(4.14)

con a = 0.948 para un nivel α = 0.05. En las Figuras 4.4 se muestra gr´aficamente los dos casos posibles. Definici´on 4.3.5. La suma acumulada de residuos recursivos cuadrados o´ CUSUMSQ se define como Pt 2 i=k Wi+1,i CU SU M SQt = PT , t ∈ [k, T ] (4.15) 2 i=k Wi+1,i

Si existe estabilidad estructural la gr´afica de CU SU M SQt versus t es una l´ınea recta. Per´ıodos en los cuales aumente la inestabilidad se notar´an apartados de esta l´ınea. Introducido por Brown et al. [1975], es una herramienta est´andar en el an´alisis de series econ´omicas.

Implementaci´on en R: la librer´ıa strucchange El objeto de esta librer´ıa son las pruebas de cambios estructurales en modelos de regresi´on. Algunas funciones de la librer´ıa

62 1. recresid: Calcula los residuales recursivos o´ errores de predicci´on a un paso estandarizados, en un modelo de regresi´on lineal. La suma acumulada de e´ stos es el CUSUM. rrc = recresid(y ˜ t + It) 2. efp: Calcula los residuales recursivos junto con la regi´on de aceptaci´on (Bandas). No es propiamente una prueba de hip´otesis sino un diagn´ostico gr´afico. Hay dos versiones. Una que utiliza los residuales recursivos y otra con base en los residuales de la regresi´on por M´ınimos Cuadrados. prueba.cusum = efp( y ˜ t + It, type = "Rec-CUSUM") plot(prueba.cusum1) prueba.cusum2 = efp(y ˜ t + It, type = "OLS-CUSUM") plot(prueba.cusum2) 3. sctest: Calcula la prueba formal de la hip´otesis nula de estabilidad estructural, y depende del tipo de residuales que se utilice para calcular. sctest(prueba.cusum1) sctest(prueba.cusum2) 4. breakpoints:Suponga el modelo Yt = X 0t β + εt , t = 1, . . . , T se asume m puntos de quiebre o de cambio en los cuales los coeficientes cambian de un valor a otro, permaneciendo constantes entre los puntos. Los puntos de quiebre son tales que Yi = X 0i βj + µi ,

ij−1 + 1 ≤ i ≤ ij ,

j = 1, . . . , m + 1.

(4.16)

Los puntos de quiebre (el ´ındice) ij se estiman con la funci´on breakpoints. bp.n = breakpoints(y2 ˜ t + It, breaks=2) summary(bp.n) # en el vector B estan las observaciones que # corresponden a cambios estructurales B= bp.n$extract.breaks(bp.n$RSS.table,breaks=2) Ejemplo 4.3.3. Con base en el Ejemplo (4.2.2), de la serie de producci´on trimestral de cemento, la gr´afica del estad´ıstico CUSUM para el modelo con tendencia lineal y variables indicadoras estacionales (5.13), pag. 86, est´a en la Figura 4.5. A partir de los resultados en la Figura 4.5 se rechazar´ıa la hip´otesis nula de estabilidad estructural en el modelo, ya que el estad´ıstico CUSUM se sale de las bandas de confianza. La prueba formal confirma el rechazo de la hip´otesis nula de estabilidad.

OLS−based CUSUM test

0.5 0.0

Empirical fluctuation process

−1.5

−4

−1.0

−0.5

0 −2

Empirical fluctuation process

1.0

1.5

Recursive CUSUM test

1960

1970

1980

1990

Time

(a) residuales recursivos

1960

1970

1980

1990

Time

(b) residuales ols

Figura 4.5: Resultados prueba CUSUM gr´afica data: prueba.cusum S = 1.703, p-value = 1.782e-05 La conclusi´on es que el modelo (5.13), pag. 86 detecta un cambio estructural en la serie. El punto de quiebre se detecta en el cuarto trimestre de 1974 (la observaci´on n´umero 76). Ver el c´odigo en R a continuaci´on para la aplicaci´on de la funci´on breakpoints. En este caso conclu´ımos que los pron´osticos con este modelo pueden no ser v´alidos. Entonces preferimos los pron´osticos con un modelo adaptativo como Holt-Winters, o´ el Modelo de Componentes Estructurales de Harvey (4.27). Tambi´en es posible ampliar el modelo de descomposici´on para inclu´ır una intervenci´on que modele el cambio en la observaci´on 76. Pero el an´alisis de intervenciones es tema de otro cap´ıtulo. El c´odigo R siguiente muestra la programaci´on de los comandos con los cuales se generaron los an´alisis y las gr´aficas de esta parte #----programa en R para analisis de estabilidad estructural #----con la serie de produccion de cemento library(strucchange) library(forecast) library(lmtest) E = read.table("F:/Curso pron´ osticos/Datos/cementq.dat", header = TRUE) attach(E)

0 −300

500

−200

−100

1000

serie

residuo

100

1500

200

300

2000

100

150

(a) Con la Serie

100 t

(b) Con los residuos

Figura 4.6: Punto de Quiebre 1974-Q4 en la Serie y en los Residuos

y = ts(y,frequency=4,start=c(1956,1),end=c(1994,3)) #----modelo con tendencia lineal y estacionalidad T = length(y) t = seq(1,T) It = seasonaldummy(y) mod1 = lm(y ˜ t + It) summary(mod1) #--- pruebe Durbin-Watson dwtest(y ˜ t + It) #----estimacion recursiva de parametros k = 2 + frequency(y) -1 + 10 n = T-k parm = mat.or.vec(n, 5) for(j in 1:n){ yj = y[1:(k+j)] tj = t[1:(k+j)] Itj = It[1:(k+j),] mod.j = lm(yj ˜ tj + Itj) parm[j,] = t(mod.j$coefficient) }

150

65 #-----grafica de las trayectorias colnames(parm)=c("beta.0", "beta.1","delta.1","delta.2", "delta.3") plot.ts(parm,main="") #-----pruebas cusum graficas prueba.cusum1 = efp(y ˜ t + It, type = "Rec-CUSUM") plot(prueba.cusum1) prueba.cusum2 = efp(y ˜ t + It, type = "OLS-CUSUM") plot(prueba.cusum2) #----pruebas cusum formales con valor p sctest(prueba.cusum1) sctest(prueba.cusum2) #----encontrar el punto de quiebre: determinar 2 puntos posibles bp.n = breakpoints(y ˜ t + It, breaks=2) summary(bp.n) B= bp.n$extract.breaks(bp.n$RSS.table,breaks=2) #---- grafica de cusum con punto de quiebre rcres = recresid(y ˜ t + It) plot(c*msum(rcres),type=’l’) abline(v=B[1],lty=2,lwd=2) abline(v=B[2],lty=2,lwd=2) r1 = mod1$residuals #----grafica de residuos OLS con punto de quiebre plot(t,r1,type=’o’,ylab=’residuo’) abline(v=B[1],lty=2,lwd=2) abline(v=B[2],lty=2,lwd=2) #----grafica de la serie con punto de quiebre

66 plot(t,y,type=’o’,ylab=’serie’) abline(v=B[1],lty=2,lwd=2) abline(v=B[2],lty=2,lwd=2) #-----------------------------------

4.4.

M´etodos de Filtrado para la Componente Estacional

En esta secci´on se exponen dos m´etodos para estimar la componente St de la serie Yt = Tt + St + εt , t = 1, . . . , T . Una con base en filtros lineales, de la forma Sˆt =

r X

uj Yt+j ,

j=−r

t = r + 1, . . ., T − r,

(4.17)

Y otra con base en un algoritmo recursivo, que suaviza la serie sin eliminar la componente estacional. Existen otras t´ecnicas de ajuste estacional que son t´ecnicas ad-hoc para calcular indices estacionales y luego usar estos para destacionalizar. Un ejemplo es el m´etodo X11-ARIMA utilizado para desestacionalizar series por el U.S. Census Bureau ahora es X12-ARIMA en www.census.gov/srd/www/x12a

4.4.1.

Algoritmo de Brockwell y Davis

Este Algoritmo est´a en en Brockwell and Davis [2002, p´ag. 24] y se describe a continuaci´on. Con base en la serie Y1 , . . . , YT , primero se estima la tendencia con una media m´ovil dise˜nada para eliminar la componente estacional y disminu´ır el ruido aleatorio, definida por: Si el per´ıodo s es par, se coloca s = 2q y para q + 1 ≤ t ≤ T − q se define 1 1 Tbt = s−1 ( Yt−q + Yt−q+1 + · · · + Yt + · · · + Yt+q−1 + Yt+q ), 2 2

(4.18)

Si el per´ıodo s es impar, se coloca s = 2q + 1, y para q + 1 ≤ t ≤ T − q, Tbt =

q X 1 Yt−j , 2q + 1

(4.19)

j=−q

Segundo, se estima Sbk para cada k = 1, 2, . . ., s. Es decir, se estima el patr´on estacional, promediando los valores de Yt − Tbt para todos los t que corresponden a la estaci´on k (son los t = sbt/sc + k). Con q = bs/2c se define

(1) Sˆk =

donde a =

(q−k)+ s

Tercero, se coloca

+1yb=

b X j=a

(T −q−k) s

Yk+sj − Tˆk+sj

k .

(4.20)

k = 1, 2, . . ., s,

(4.21)

b−a+2

1 X ˆ(1) (1) Sˆk = Sˆk − Sk , q k=1

P para garantizar sk=1 Sˆk = 0 que es una condici´on equivalente a eliminar la u´ ltima estaci´on, y tiene como objeto garantizar la identificabilidad de la componente St . Finalmente se repite el patr´on obtenido el n´umero de veces necesario para formar Sˆt, t = 1, . . . , T. Ejemplo 4.4.1. Tomando por ejemplo T = 250, s = 12 y q = 6 entonces

1 Tˆt = 12

1 1 Yt−6 + Yt−5 + · · · + Yt+5 + Yt+6 , 2 2

con 7 ≤ t ≤ 244, entonces se tiene que

(1) Sˆk =

b X j=a

Yk+12j − Tˆk+12j b−a+2

donde

(6 − k)+ a= + 1, 12 (250 − 6 − k) b= . 12 Para k = 1 20

1 X (1) Sˆ1 = Y1+12j − Tˆ1+12j 21 j=1 1 = (Y1+12 − Tˆ1+12 ) + (Y1+24 − Tˆ1+24 ) + · · · + (Y1+240 − Tˆ1+240) . 21

68 Note que Tˆ244 es el u´ ltimo valor de Tˆt, usamos Tˆ241 con Y241 . Para k = 12 19 19 1 X 1 X (1) ˆ ˆ Y12+12j − T12+12j = Y12(1+j) − Tˆ12(1+j) S12 = 20 20 j=0 j=0 1 = (Y12 − Tˆ12) + (Y24 − Tˆ24 ) + · · · + (Y240 − Tˆ240 ) . 20

El algoritmo se puede implementar en una funci´on en R como sigue. seascomp = function(x,p){ #----------algoritmo Brockwell-Davis n 1, 2 2 4 luego Yt es invertible.

6.2.3.

Funci´on fac parcial de un Proceso MA(q) invertible

Suponga un proceso Yt ∼ M A(q) invertible, Yt = θq (L)(εt).

(6.7)

Considere θq (z) = 1 + θ1 z + · · · + θq z q entonces θq (z) 6= 0, |z| ≤ 1, luego la funci´on tiene desarrollo es serie de Taylor alrededor de z = 0, dado por

1 θq (z)

∞

X 1 = 1 + ψ1 z + ψ2 z 2 + . . . = ψj z j , θq (z)

ψ0 = 1,

(6.8)

j=0

con

P∞

1 θq (L)

ψ 2 < ∞, donde ψj → 0 si j → ∞. Multiplicando ambos miembros de (6.7) por se obtiene j=0

εt =

1 Yt = ψ(L)(Yt) = Yt + ψ1 Yt−1 + ψ2 Yt−2 + . . . θq (L)

(6.9)

Y despejando Yt Yt = −ψ1 Yt−1 − ψ2 Yt−2 − · · · + εt ,

(6.10)

de donde conclu´ımos que si hace la regresi´on de Yt sobre los primeros k rezagos Yt−j , j = 1, . . . , k, entonces el k-´esimo coeficiente es α(k) = ψ(k) 6= 0, ∀k y como ψ(k) → 0 entonces α(k) → 0 cuando k → ∞. Por tanto, la facp de un MA(q) decrece a cero. En las Figuras siguientes 6.3 se observa la fac y la facp de un MA(3).

True PACF

0.6 −0.2 0.2

True PACF

0.8 0.4 0.0

True ACF

1.0

True ACF

Lag

(a) F AC

Lag

(b) F ACP

Figura 6.3: F AC y F ACP de un M A(3).

6.2.4.

Implementaci´on en R

En R para identificar se usan las funciones acf y pacf y para estimar una de las funciones usadas es arma de la librer´ıa tseries.

Ejemplo 6.2.3. library(forecast,tseries) n = 300 theta = c(-1,-0.4,-0.4) (Mod(polyroot(theta))) y = arima.sim(list(order=c(0,0,2), ma=theta[2:3]), n=n, sd=sqrt(2.3)) layout(1:3) ts.plot(y) acf(y,30) pacf(y,30) # Estimaci´ on: Funci´ on arma librer´ ıa tseries modelo.y = arma(x=y, order=c(0,2)) summary(modelo.y)

100

0.2 0.4 −0.1

0.6

Partial ACF

Series y

0.0

ACF

Series y

Lag

(a) F AC

(b) F ACP

Figura 6.4: F AC y F ACP del Ejemplo. pred.y = predict(modelo.y, n.ahead=2) plot(seq(1,9,1), c(tail(y), pred.y$pred), type=’b’) points(seq(7,9,1), pred.y$pred, type=’b’, col=’red’)

6.3.

Procesos Autorregresivos de Orden p, AR(p)

Definici´on 6.3.1 (Proceso AR(p)). Se dice que Yn , n ∈ Z sigue un proceso AR(p) de media cero si Yn = ϕ1 Yn−1 + ϕ2 Yn−2 + · · · + ϕp Yn−p + εn , (6.11) donde εn ∼ RB(0, σ 2 ) y p = 1, 2, . . .. Usando el operador de rezago L se puede escribir (6.11) como ϕp(L)(Yn ) = εn , (6.12) con ϕp (z) = 1 − ϕ1 z − ϕ2 z 2 − · · · − ϕpz p , z ∈ C, el polinomio autorregresivo.

6.3.1.

Condici´on Suficiente para que un AR(p) sea Estacionario en Covarianza

La condici´on suficiente para que Yt ∼ AR(p) sea estacionario en covarianza es que las p ra´ıces de la ecuaci´on ϕp(z) = 0, zi , i = 1, 2 . . ., p cumplan |zi | > 1.

(6.13)

q N´otese que si zj = aj ± ibj entonces |zj | = a2j + b2j . La condici´on (6.13) no es, sin embargo, necesaria. En palabras, la condici´on (6.13) se describe como “ para que un proceso

101 autorregresivo de orden p sea estacionario en covarianza, es suficiente que las ra´ıces del polinomio autorregresivo est´en por fuera del c´ırculo unitario”. El c´ırculo unitario aparece en la Figura 6.5. En esta figura se observa la posici´on de la ra´ız zj y su conjugado z¯j .

Figura 6.5: C´ırculo Unitario Si el proceso Yt es estacionario en covarianza se cumple que su media es constante, Yt ≡ µ, y tiene funci´on de autocovarianza, Cov(Yt , Yt+k ) = R(k), k ∈ Z. De ahora en adelante siempre se asumir´a que los procesos AR(p) cumplen la condici´on (6.13).

6.3.2.

Algunas Propiedades de los Procesos AR(p) Estacionarios

Proposici´on 6.3.1. Para un proceso Yt ∼ AR(p), (6.11) se cumple: i) E(Yt) = 0 y ii) Pp j=1 ϕj < 1.

Demostraci´on. i). Si Yt es estacionario en covarianza entonces E(Yt ) = µ. Adem´as, E(Yt) = ϕ1 E(Yt−1 ) + ϕ2 E(Yt−2 ) + · · · + ϕp E(Yt−p) + 0, pero todas las esperanzas son µ luego µ = ϕ1 µ + ϕ2 µ + · · · + ϕpµ.

Si µ 6= 0 entonces 1 = ϕ1 + · · · + ϕp , por tanto el polinomio autorregresivo evaluado en z = 1 se anula, es decir, ϕp(1) = 0, lo cual es una contradicci´on ya que ∀ z ∈ C, |z| ≤ 1 ϕp(z) 6= 0. Luego debe tenerse que µ = 0, es decir, el proceso definido en (6.11) es de Pp media cero. Para ii), si se asume que j=1 ϕj ≥ 1, el caso = 1 implica nuevamente que P ϕp(1) = 0, lo cual es una contradicci´on. Si se asume pj=1 ϕj > 1 entonces ϕp (1) < 0. Pero ϕp(0) = 1. Como la funci´on ϕp(x) es continua en [0, 1] debe existir un punto 0 < x < 1 en

102 el cual ϕp (x) = 0. Pero eso es nuevamente una contradicci´on porque este x es una ra´ız de ϕp (z) = 0 que cumple |z| < 1, lo cual no puede ocurrir. Un proceso Yt ∼ AR(p) con E(Yt ) = µ 6= 0 se define mediante la relaci´on siguiente. ϕp (L)(Yt) = ϕ0 + εt ,

(6.14)

donde ϕ0 = ϕp (L)(µ) = (1 − ϕ1 − ϕ2 − · · · − ϕp )µ. N´otese que tambi´en se puede escribir Yt = (1−ϕ1 −· · ·−ϕp )µ+ϕ1 Yt−1 +· · ·+ϕp Yt−p +εt , de donde Yt − µ = ϕ1 (Yt−1 − µ) + · · · + ϕp (Yt−p − µ) + εt . Es decir, el proceso (Yt − µ) es AR(p) de media cero.

La Funci´on de Autocovarianza de los Procesos AR(p) La funci´on de autocovarianza de un proceso Yt ∼ AR(p), R(k) = Cov(Yt , Yt+k ), se puede calcular resolviendo una ecuaci´on recursiva lineal denominada la ecuaci´on de Yule–Walker. P Proposici´on 6.3.2. Suponga un proceso AR(p), Yn = pj=1 ϕj Yn−j + εt , que satisface la condici´on ((6.13)). Su funci´on de autocovarianza R(k) satisface la ecuaci´on recursiva R(k) =

p X j=1

ϕj R(k − j), k = 1, 2, . . . .

(6.15)

denominada ecuaci´on de Yule–Walker. P Demostraci´on. Colocando µ = E(Yn ), como Yn = pj=1 ϕj Yn−j + εn , al tomar esperanza P en ambos miembros de esta identidad se obtiene µ = pj=1 ϕj µ+0. Restando las expresiones P anteriores se obtiene Yn − µ = pj=1 ϕj (Yn−j − µ) + εn . Multiplicando ambos miembros de esta identidad por Yn−k − µ, con k ≤ n, y tomando valor esperado E(.) se obtiene R(k) = E((Yn − µ)(Yn−k − µ)) p X = ϕE((Yn−j − µ)(Yn−k − µ)) + E(εn (Yn−k − µ)) j=1 p

X j=1

ϕj R(k − j).

103 En el resultado anterior se cumple E(εn (Yn−k − µ)) = 0 porque, a partir de la definici´on del proceso Yn en (6.11), Yn−k depende de εs con s ≤ n−k, que son variables incorrelacionadas con εn .

La Varianza de los Procesos AR(p) Si Yt ∼ AR(p) de media cero, ϕp(L)(Yn ) = εn , para εn ∼ RB(0, σ 2 ). Adem´as, se asumi´o que se cumple que ∀z, |z| ≤ 1, ϕp(z) 6= 0. Entonces el cociente ϕp1(z) se puede desarrollar en serie de potencias de z, y colocar ∞

X 1 = ψj z j , ϕp (z) j=0

para ciertos coeficientes (ψj , j = 0, 1, . . .), con ψ0 = 1. Por tanto, se puede colocar Yn =

1 (εn ) = εn + ψ1 εn−1 + ψ2 εn−2 + . . . . ϕp(L)

(6.16)

P 2 Tomando varianza a ambos miembros de (6.16), se obtiene V ar(Yn ) = σ 2 ∞ j=0 ψj . La identidad (6.16) se puede interpretar como que todo proceso AR(p) que cumple la condici´on (6.13) es equivalente a un M A(∞). Tambi´en se denomina (6.16) la representaci´on causal de Yt . Los coeficientes ψj se pueden calcular con la funci´on en R, ARMAtoMA(ar = numeric(0), ma = numeric(0), lag.max). En el argumento “ar=” se ingresa el vector de coeficientes autorregresivos (ϕ1 , ϕ2 , . . . , ϕp). El argumento “lag.max” es el n´umero de coeficientes ψj deseados.

Estimaci´on de la FAC de un AR(p) La fac de un AR(p), ρ(k) = Corr(Yt , Yt+k ), siguiente.

k = 1, 2, . . . , p, p + 1, . . . , cumple lo

1. Dadas la matriz A y los vectores ϕ, ρ 

   A=   

1 ρ(1) ρ(2) .. .

ρ(1) ρ(2) ρ(3) .. .

··· ··· ···

ρ(p − 1) ρ(p − 2) · · ·

ρ(p − 1) ρ(p − 2) ρ(p − 3) .. . 1



      , ϕ =      

ϕ1 ϕ2 .. . ϕp





    , ρ =     

ρ(1) ρ(2) .. . ρ(p)



  ,  

104 entonces es v´alido el sistema lineal p × p Aϕ = ρ.

(6.17)

Luego, dados los estimadores de la fac ρˆ(1), . . . , ρˆ(p) se puede resolver (6.17) y colocar ϕ ˆ = Aˆ−1 ρˆ. Estos estad´ısticos ϕ ˆ as´ı calculados se denominan los estimadores de Yule-Walker de ϕ 2. Se puede comprobar que la fac de un AR(p) cumple la ecuaci´on siguiente. ρ(k) = ϕ1 ρ(k − 1) + ϕ2 ρ(k − 2) + · · · + ϕp ρ(k − p),

k = p, p + 1, . . . (6.18)

Entoces (6.18) forma una ecuaci´on en diferencias finitas con condiciones iniciales ρ(1), . . ., ρ(p), para ρ(k), k ≥ p + 1, con soluci´on ρ(k) = s1 g1k + s2 g22 + · · · + sp g2p,

(6.19)

donde gi = 1/zi y zi es la i-´esima ra´ız de la ecuaci´on caracter´ıstica 1 − ϕ1 z − ϕ2 z 2 − · · · − ϕp z p = 0

(6.20)

con |zi | > 1 ⇔ |gi | < 1, luego se debe cumplir que ρ(k) → 0, k → ∞ Nota 6.3.1. Si gi ≈ 1, por ejemplo gi = 1 − ε se tendr´a si gik = si (1 − ε)k y por tanto ρ(k) decae a cero m´as lento que si gi = ε.

0.8 0.4 0.0

0.4

True ACF

0.8

True ACF

0.0

True ACF

Lag

10 Lag

(a) ϕ = 0.35

(b) ϕ = 0.88

Figura 6.6: FAC de Yt = ϕYt−1 + εt .

105 FACP de los Procesos AR(p) La FACP de un procesos AR(p) es α(k) tal que α(k) ˆ es el coeficiente βˆk,k en la regresi´on Yt = β0 + βk,1 Yt−1 + · · · + βk,k Yt−k + at ,

k=2

(6.21)

pero como βk,k = 0 si k ≥ p + 1 entoces α ˆ (k) = 0 si k ≥ p + 1

0.8 0.4 0.0

True ACF

Lag

Figura 6.7: FACP Muestral de AR(p). Ejemplo 6.3.1. Sea Yt ∼ AR(2) con Yt = ϕ1 Yt−1 + ϕ2 Yt−2 + εt , Yt = 1.5Yt−1 − 0.9Yt−2 + εt ,

i.i.d.

εt ∼ N (0, σ 2)

ϕ2 (z) = 1 − 1.5z + 0.9z 2 = 0 z = 0.83 ± 0, 64i,

ecuaci´on caracter´ıstica

|z| = 1.054 > 1

luego Yt es estacionario en covarianza, adem´as ρ(k) = 1 − 1.5ρ(k − 1) + 0.9ρ(k − 2), k ≥ 2 ϕ1 1.5 ρ(0) = 1, ρ(1) = = = 0.789. 1 − ϕ2 1.9

6.3.3.

Proceso AR(1)

El proceso AR(1) se ha utilizado anteriormente por ejemplo, en la prueba Durbin-Watson. A continuaci´on se desarrollan algunas de sus propiedades. Si Yt es un AR(1) de media µ, entonces est´a definido por Yt = ϕ0 + ϕ1 Yt−1 + εt ,

ϕ0 = (1 − ϕ1 )µ,

(6.22)

106 donde el proceso Yt es estacionario si todas la ra´ıces z de la ecuaci´on ϕ1 (z) = 0 caen fuera del circulo unitario |z| > 1. Luego el AR(1) es estacionario en covarianza si y solo si |ϕ1 | < 1. Definici´on 6.3.2 (Marcha Aleatoria). Se dice que Yt es una marcha aleatoria (Random Walk) si cumple Yt = µ + Yt−1 + εt , (6.23) n´otese que es un AR(1) con ϕ1 = 1.

Propiedades del AR(1) 1. E(Yt) = µ =

ϕ0 1 − ϕ1

2. Cov(Yt , Yt+k ) = 3. ρ(k) = ϕk1 ,

σ 2 ϕk1 , 1 − ϕ2 1

k = 0, 1, . . .

−1 < ϕ1 < 1.

Nota 6.3.2. Diebold [1999, p´ag. 138], Si Yt ∼ AR(1) de media cero estacionario en covarianza entonces Yt = ϕ1 Yt−1 + εt ,

εt ∼ R.B.(0, σ 2)

es decir (1 − ϕ1 L)Yt = εt y se puede escribir como Yt =

1 εt , 1 − ϕ1 L

si se cumple que ∞

f (z) =

X j 1 = ϕ1 z j = 1 + ϕ1 z + ϕ21 z 2 + . . . 1 − ϕ1 z j=0

por que |ϕ1 z| < 1 ya que |ϕ1 | < 1 y |z| < 1. Entonces Yt = εt + ϕ1 εt−1 + ϕ21 εt−2 + . . . ,

107 y como los εt son incorrelacionados V ar(Yt ) = σ 2 + ϕ1 σ 2 + ϕ21 σ 2 + . . . 2

= σ (1 + ϕ1 +

ϕ21

+ . . .) = σ

1 1 − ϕ1

Varianza Incondicional

Nota 6.3.3. Si Yt ∼ AR(1) de media cero estacionario en covarianza entonces E(Yt|Yt−1 ) = E(ϕ1Yt−1 + εy |Yt−1 ) = ϕ1 Yt−1 , y si se asume que εt son independientes de yt−1 , Yt−2 , . . . V ar(Yt |Yt−1 ) = V ar(ϕ1 Yt−1 + εt |Yt−1 )

= ϕ21 V ar(Yt−1|Yt−1 ) + V ar(εt) = σ 2

Varianza Condicional

6.4. Procesos Autoregresivos y de Medias M´oviles ARMA(p,q) Si en un proceso Yt ∼ AR(p) Yt − ϕ1 Yt−1 − ϕ2 Yt−2 − · · · − ϕp Yt−p = εt ,

εt ∼ R.B.(0, σ 2),

(6.24)

se cambia εt por un modelo M A(q), Zt = εt + θ1 εt−1 + · · · + θq εt−q entonces (6.24) queda Yt − ϕ1 Yt−1 − · · · − ϕp Yt−p = εt + θ1 εt−1 + · · · + θq εt−q ,

(6.25)

donde εt ∼ RB(0, σ 2 ). El efecto de este cambio es que los errores no se toman incorrelacionados sino con autocorrelaci´on d´ebil. Se define entonces un proceso ARMA(p,q) como un modelo que combina las propiedades de memoria larga de los AR(p) con las propiedades de ruido d´ebilmente autocorrelacionado en los MA(q), y que tiene suficiente flexibilidad y parsimonia. Usando la notaci´on del operador de rezago (6.25) se puede definir el proceso ARM A(p, q) por Definici´on 6.4.1. Un proceso Yt ∼ ARM A(p, q) se define mediante la ecuaci´on (6.25), o tambi´en por ϕp (L)(Yt) = θq (L)(εt),

t ∈ Z,

(6.26)

P P donde εt ∼ RB(0, σ 2), y ϕp (z) = 1 − pj=1 ϕj z j , θq (z) = 1 + qj=1 θj z j son los polinomios autoregresivo y de media m´ovil respectivamente.

108 Las condiciones de estacionariedad de la parte AR(p) y de invertibilidad de la parte MA(q) se asumen en el modelo ARMA(p,q), (6.26). Por lo tanto, se asume que las ra´ıces de las ecuaciones ϕp (z) = 0 y θq (z) = 0 est´an fuera del c´ırculo unitario. Adem´as se asume que estos polinomios no tienen ra´ıces en com´un. Si se cumplen estas condiciones el proceso Yt ∼ ARM A(p, q) es estacionario e identificable. Ejemplo 6.4.1. Sea Yt ∼ ARM A(1, 1) dado por ϕ1 (L)(Yt) = θ1 (L)(εt)

(6.27)

donde ϕ1 (L) = 1 − ϕL y θ1 (L) = 1 + θL. Es decir Yt = ϕYt−1 + εt + θεt−1 . Si |ϕ| < 1 y |θ| < 1 es estacionario e invertible. Por ejemplo Yt = 0.9Yt−1 + εt − 0.4εt−1 , con εt ∼ RB(0, σ 2) Ejemplo 6.4.2. Consideremos un modelo de descomposici´on con tendencia lineal y estacionalidad de per´ıodo 12, modelada por variables indicadoras donde el residuo estructural es un proceso ARM A(1, 1). Entonces el modelo se escribe como el sistema de dos ecuaciones siguiente. Yt = β0 + β1 t +

11 X

δj Ij (t) + εt ,

j=1

εt = ϕεt−1 + at + θat−1 , con at ∼ RB(0, σ 2 ), el residuo del proceso arma. N´otese que el n´umero de par´ametros de este modelo es 16, incluyendo la varianza σ 2 . Ejemplo 6.4.3. Considere el proceso Yt ∼ ARM A(2, 1) dado por Yt = 2 +

1 − 0.4L εt , 1 − 1.5L + 0.9L2

εt ∼ R.B.(0, σ 2)

osea Yt = 2(1 − 1.5 + 0.9) + 1.5Yt−1 − 0.9Yt−3 + εt − 0.4εt−1 con ecuaci´on caracter´ıstica 1 − 1.5z + 0.9z 2 = 0

109 y sus ra´ıces dadas por z=

1.5 ±

1.52 − 4(0.9) = 0.83 ± 0.645i 2(0.9)

por tanto |z| = 1.05 > 1 es un proceso estacionario en covarianza e invertible.

6.4.1.

Propiedades de los Modelos ARMA

1. Suponga Yt ∼ ARM A(p, q) entonces E(Yt) = 0. Si el proceso es estacionario P P∞ j entonces se puede expresar θq (z)/ϕp(z) = ∞ j=0 ψj z con ψ0 = 1 y j=0 |ψj | < ∞. A partir de ϕp(L)Yt = θq (L)εt, se puede escribir ∞

Yt =

X θq (L) εt = ψj εt−j = εt + ψ1 εt−1 + ψ2 εt−2 + . . . , ϕp (L) j=0

por tanto E(Yt) = E(εt + ψ1 εt−1 + ψ2 εt−2 + . . . ) = 0. 2. En caso de ser E(Yt) = µ 6= 0 se coloca Yt = µ +

θq (L) εt ϕp (L)

(6.28)

de donde ϕp(L)Yt = ϕp(1)µ + θq (L)εt. Por ejemplo, sea Yt ∼ ARM A(1, 1) con E(Yt) = µ entonces Yt = µ +

1 + θL εt 1 − ϕL

luego (1 − ϕL)Yt = (1 − ϕL)µ + (1 − θL)εt pero (1 − ϕL)µ = µ − ϕµ = (1 − ϕ)µ, luego Yt = (1 − ϕ)µ + ϕYt−1 + εt + θεt−1 .

110 3. La funci´on de autocovarianza de un proceso Yt ∼ ARMA(p,q) estacionario de media cero. Si se indica por R(k) = Cov(Yt , Yt+k ) su funci´on de autocovarianza, para k = 0, 1, . . . un m´etodo para calcular esta funci´on se basa en la representaci´on P j ϕp(L)Yt = θq (L)εt, con θq (z)/ϕp(z) = ∞ j=0 ψj z . Multiplicando ambos miembros por Yt−k y tomando esperanza E(.) se obtienen las ecuaciones recursivas siguientes, similares a las ecuaciones Yule-Walker para AR(p), (6.15). Defina n = max(p, q + 1),  σ 2 Pq θ ψ , si k = 0, 1, . . ., n − 1 j=k j j−k R(k)−ϕ1 R(k−1)−. . .−ϕp R(k−p) = 0, si k = n, n + 1, . . .. (6.29) Ejemplo 6.4.4. (tomado de Brockwell and Davis [2002], pag. 93). Considere el proceso ARMA(2,1) dado por (1−L+ 14 L2 )Yt = (1+L)εt. Entonces n = max(p, q+ 1) = 2, por tanto, para k = 0, 1 se tiene el sistema lineal

1 R(0) − R(1) + R(2) = σ 2 (ψ0 θ0 + ψ1 θ1 ) = σ 2 (1 + ψ1 ), 4 1 R(1) − R(0) + R(1) = σ 2 (θ1 ψ0 ) = σ 2 ψ0 = σ 2 . 4

(6.30)

Para calcular los coeficientes (ψj , j = 0, 1, . . .) se puede utilizar la funci´on de R, ARMAtoMA(a,b,m), la cual calcula los coeficientes ψj , j = 1, 2, . . ., m, dados los vectores a = (ϕ1 , . . ., ϕp) y b = (θ1 , . . . , θq ). Entonces, escribiendo ARMAtoMA(c(1.0, -0.25), 1.0, 10), se obtiene el vector [1] 2.00000000 1.75000000 1.25000000 0.81250000 0.50000000 [6] 0.29687500 0.17187500 0.09765625 0.05468750 0.03027344

de donde ψ1 = 2. Usando la segunda ecuaci´on de (6.29), se obtiene R(2) = R(1) − 1 4 R(0), luego, reemplazando en el sistema (6.30), y resolviendo, se obtienen R(0) = 32σ 2 /3, R(1) = 28σ 2 /3. Utilizando la segunda ecuaci´on de (6.29), 1 R(k) = R(k − 1) − R(k − 2), k = 2, 3, . . . 4 se puede calcular la autocovarianza recursivamente. 4. La varianza de un proceso Yt ∼ ARMA(p,q) estacionario de media cero. Es evidente que a partir de la primera ecuaci´on en (6.29) se puede definir un sistema lineal una de cuyas inc´ognitas es R(0), la cual se resuelve en funci´on de σ 2 .

111

6.4.2.

Librer´ıas para identificaci´on, estimaci´on y pron´osticos de modelos ARMA

El plan de an´alisis con procesos ARMA consiste en

1. Identificar el modelo ARMA(p,q). 2. Estimar el modelo. 3. Chequeo del ajuste del modelo a los datos. 4. Pron´osticos con el modelo. O simulaci´on del modelo.

Algunas de las librer´ıas y funciones a utilizar para este an´alisis son

1. stat, con la funci´on arima(), estima primero por m´ınimos cuadrados condicionales y luego por m´axima verosimilitud. 2. tseries, con la funci´on arma(), estima mediante m´ınimos cuadrados condicionales. 3. forecast, con la funci´on auto.arima(), para identificaci´on de modelos ARIMA. 4. FitAR, FitARMA, para estimaci´on de AR(p) y ARMA(p,q). 5. timsac, con la funci´on autoarmafit(), para identificaci´on del modelo ARMA(p,q) con menor AIC. El programa asume que el modelo ARMA(p,q) se define con θq (z) = P 1 − qj=1 θj z j , es decir, los coeficientes θj se cambian de signo en esta librer´ıa. Tambi´en tiene la funci´on armafit() para ajuste de modelos ARMA.

6.4.3.

Identificaci´on de modelos ARMA

No es f´acil identificar los modelos ARM A(p, q) mediante la FAC y FACP. Si (q ≥ p ≥ 1) entonces para (k ≥ q + 1) se cumple (6.18) y para 1 ≤ k ≤ p − 1, ρ(k) no tiene un patr´on general, luego la FAC muestral presenta un patr´on definido solamente para k ≥ q + 1.

112

Figura 6.8: Fac Muestral de ARM A(p, q).

Una alternativa consiste en buscar una pareja de o´ rdenes (p, q) dentro de un rango inicial, por ejemplo p, q = 0, 1, . . ., 10, que minimice alguna funci´on de εˆ2 , como por ejemplo el AIC o´ el BIC.

6.4.4.

Estimaci´on de modelos ARMA

La estimaci´on de procesos Yt ∼ ARMA(p,q) se basa en un supuesto: que el vector Y = (Y1 , . . . , Yn )0 se distribuye Normal multivariado con media µ, y matriz de covarianzas Σ = [Cov(Yi , Yj )]n×n . Como el proceso se asume estacionario, se cumple Cov(Yi , Yj ) = R(j − i), donde R(k) es la funci´on de autocovarianza de Yt . La forma de Σ es la de una matriz tipo Toeplitz: las diagonales descendentes son constantes: 

   Σ=   

R(0) R(1) R(2) .. .

R(1) R(0) R(1) .. .

··· ··· ···

R(n − 1) R(n − 2) · · ·

R(n − 1) R(n − 2) R(n − 3) .. . R(0)



   .   

(6.31)

Por ejemplo, para Yt un AR(p), R(k) se calcula mediante las ecuaciones Yule-Walker, (6.15), P R(k) = µ + pj=1 ϕj R(k − j). Por tanto, colocando β = (µ, σ 2 , ϕ1 , . . . , ϕp, θ1 , . . ., θq )0 , la matriz Σ depende del vector β, y se escribe Σ(β). Este supuesto permite implementar la estimaci´on por m´axima verosimilitud. Se escribe la densidad Normal Multivariada como 1 1 −1 0 q f (y, β) = exp − (y − µ)Σ(β) (y − µ) 2 (2π)n/2 det(Σ(β))

(6.32)

113 donde µ = (µ, . . . , µ)0 ∈ Rn . La funci´on de log-verosimilitud se define a partir del logaritmo de la densidad (6.32), log(f (y, β)), y est´a dada por L(β) :=

n X

log f (yj , β)

j=1

1 n 1 − log(2π) − log det(Σ(β)) − (y − µ)0 Σ(β)−1 (y − µ). 2 2 2

(6.33)

El estimador ML de m´axima verosimilitud de β se define como b = argmin(−L(β)). β

(6.34)

La identificaci´on con base en el AIC se basa en comparar varios modelos ARM A(p, q) para valores de, por ejemplo, p = 0, 1, . . ., 20 y p = 0, 1, . . ., 20 con base en el criterio AIC el cual est´a definido por b + 2k, AIC = −2L(β)

k = p+q

(6.35)

Entonces identifica un posible modelo parsimonioso escogiendo el modelo con el m´ınimo AIC. Ejemplo 6.4.5. Suponga que se requiere simular un ARM A(3, 2) tal que se escogen las ra´ıces de ϕ3 (z) y θ2 (z), con m´odulos fuera del c´ırculo unitario. #-------------- simulacion de un ARMA(3,2) #-------------- defina las raices del polinomio autorregresivo grado 3 z1 = complex(3) z1[1] = -0.8 - 1.3i z1[2] = conj(z1[1]) z1[3] = 1.2 Entonces las tres ra´ıces son z1[1], z1[2], z1[3] y los coeficientes del polinomio autorregresivo de grado 3, m´onico, φ3 (z), φ3 (z) = (z − z1[1])(z − z1[2])(z − z1[3]) = 1 − ϕ1 z − ϕ2 z 2 − ϕ3 z 3

= 1 − 0.1466381 ∗ z − 0.1430615 ∗ z 2 − 0.3576538 ∗ z 3 se calculan a continuaci´on en el vector a.

114 a = poly.calc(z1) a = a/a[1] #-------------- defina las raices del polinomio de media movil de grado 2 z2 = complex(2) z2[1] = -1.2 -0.5i z2[2] = Conj(z2[1]) #-------------- los coeficientes estan en el vector b b = poly.calc(z2) (b = b/b[1]) #------------- usar la funcion arima.sim con los coeficientes a y b. n = 300 y = arima.sim(list(order=c(3,0,2), ar=-a[2:4], ma=b[2:3]), n=n, sd=sqrt(0.3)) #-------------- verificar el modelo require(forecast) auto.arima(y) #-------------- auto.arima identifica un arma(3,2) #-------------- para su estimacion se usa la instruccion mod1 = arima(y1, c(3,0,2)) #-------------- los pronosticos con este arma(2,3) py1 = predict(mod1, n.ahead=20) plot(1:70, c(y[(3200-49):3200],py1$pred), ylim=c(min(py1$pred-1.64*py1$se),max()), type=’b’, col=2) points(51:70,py1$pred, type=’b’, col=’blue’) points(51:70, py1$pred+1.64*py1$se, type=’l’, col=’blue’) points(51:70, py1$pred-1.64*py1$se, type=’l’, col=’blue’)

6.5.

Modelos ARMA Estacionales, SARMA

Se define un modelo ARMA Estacional (1 ) como un modelo ARMA que incluye un modelo ARMA para los tiempos dentro de cada per´ıodo estacional. Y otro modelo ARMA entre las estaciones, por ejemplo, entre los eneros, febreros, etc.. Se escribe Yt ∼ 1

http://robjhyndman.com/researchtips/cyclicts/

115 ARM A(p, q)(ps, qs )[s] . Yt ∼ ARM A(p, q)(ps, qs )[s] ⇔ Φp(L)Φps (L)Yt = Θq (L)Θqs (L)t.

(6.36)

donde (p, q) son los o´ rdenes AR y MA del modelo ARMA ordinario y (ps , qs) los o´ rdenes AR y MA estacionales, entre estaciones. Por ejemplo, el modelo Yt ∼ ARM A(1, 1)(1, 2)[s] se indica por (1 − ϕ1 L)(1 − ϕs Ls )Yt = (1 + θ1 L)(1 + θs Ls + θ2s L2s )t ,

(6.37)

donde, por ejemplo, s = 4, 12. N´otese que (6.37) es un tipo particular de ARMA. Otro ejemplo es Yt ∼ ARM A(1, 2)(1, 1)[s] que se indica por (1 − ϕ1 L)(1 − ϕs Ls )Yt = (1 + θ1 L + θ2 L2 )(1 + θs Ls )t,

(6.38)

Una clave para la identificaci´on de modelos SARMA es que si se observa por ejemplo (6.38), con s = 12, es un ARMA en el cual hay un coeficiente autorregresivo ϕ1 ϕ12 en el rezago 1 + 12 = 13, y tambi´en hay un coeficiente de media m´ovil θ1 θ12 en el rezago 1 + 12 = 13. Esto podr´ıa verse en la fac y en la fac parcial. La librer´ıa library(TSA) tiene la funci´on armasubsets que calcula un panel para identificaci´on para

Identificaci´on y Estimaci´on de los SARMA Los pasos siguientes estiman el modelo de tendencia lineal inicial. Luego muestra un procedimiento de identificaci´on gr´afico para el SARMA. De lo que se trata es de escoger los o´ rdenes p, q, ps, qs . Por ejemplo, si se supone que la serie tiene tendencia lineal y componente estacional de per´ıodo s = 12. Para identificar un posible modelo SARMA en los residuos se programan las instrucciones siguientes. #---------------------mod2 = lm(yi ˜ t ) summary(mod2) r = mod2$residuals # identificador grafico library(TSA) res=armasubsets(y=r,nar=14,nma=14,y.name=’r’, ar.method=’ols’) plot(res)

116

−1100

La u´ ltima instrucci´on grafica el panel (2,2) de la Figura 6.9. En la Figura 6.9 aparecen las

−7500 −7500 −7500 −7500 −7500

−7500 −7500 −7000 −5400

−4700 −4700 −4700 −4700 −4700

−4700 −4700 −4600

−1100

−1100 −1100

−590

−1000

−8500 −8400 −8300 −8100 −7800 −7300 −6200 −5400

En la celda (1,2), arriba-derecha, se observan rezagos autorregresivos p = 5,7, y de media m´ovil q = 2,4. Entonces un modelo factible es un ARMA(4,4). N´otese que con auto.arima() se identifica un ARMA(4,4).

En la gr´afica (1,1) aparecen muy evidentes los rezagos autorregresivos hasta p=3. Entonces se puede identificar un AR(3).

En cada gr´afica la informaci´on sobre los rezagos a inclu´ır en el modelo est´a orgarnizada en dos sectores. A la izquierda los autorregresivos, identificada con “r-lag” y a la derecha los de media m´ovil, “error-lag”. Los modelos m´as opcionados o´ factibles son los que incluyen los rezagos con menor BIC, es decir, los cuadros oscuros de la primera l´ınea.

identificaciones de un AR en la celda (1,1). En las celdas (1,2) y (2,1) aparecen modelos ARMA. Y en la (2,2) un SARMA.

Figura 6.9: Posibles submodelos (1,1)=AR, (1,2),(2,1)=ARMA, (2,2)=SARMA

−3600

(Intercept) r−lag1 r−lag2 r−lag3 r−lag4 r−lag5 r−lag6 r−lag7 r−lag8 r−lag9 r−lag10 r−lag11 r−lag12 r−lag13 r−lag14 error−lag12 error−lag13 error−lag14 error−lag1 error−lag2 error−lag3 error−lag4 error−lag5 error−lag6 error−lag7 error−lag8 error−lag9 error−lag10 error−lag11 (Intercept) r−lag1 r−lag2 r−lag3 r−lag4 r−lag5 r−lag6 r−lag7 r−lag8 r−lag9 r−lag10 r−lag11 r−lag12 r−lag13 r−lag14 error−lag4 error−lag5 error−lag6 error−lag7 error−lag8 error−lag9 error−lag10 error−lag11 error−lag12 error−lag13 error−lag14 error−lag1 error−lag2 error−lag3

(Intercept) r−lag1 r−lag2 r−lag3 r−lag4 r−lag5 r−lag6 r−lag7 r−lag8 r−lag9 r−lag10 r−lag11 r−lag12 r−lag13 r−lag14 error−lag2 error−lag3 error−lag4 error−lag5 error−lag6 error−lag7 error−lag8 error−lag9 error−lag10 error−lag11 error−lag12 error−lag13 error−lag14 error−lag1 (Intercept) r−lag1 r−lag2 r−lag3 r−lag4 r−lag5 r−lag6 r−lag7 r−lag8 r−lag9 r−lag10 r−lag11 r−lag12 r−lag13 r−lag14 error−lag1 error−lag2 error−lag3 error−lag4 error−lag5 error−lag6 error−lag7 error−lag8 error−lag9 error−lag10 error−lag11 error−lag12 error−lag13 error−lag14

BIC BIC

117 En la celda (2,1), abajo-izquierda, se observan rezagos autorregresivos p = 4. Y se observan rezagos de media m´ovil q = 4,7. Un posible modelo es un ARMA(4,4). Con auto.arima() se identifica un ARMA(4,5). En la celda (2,2), abajo-derecha, se observa algo que no aparece en las anteriores. Aparecen rezagos autorregresivos 1,12,13. Lo que d´a lugar a suponer una estructura (1 − ϕL)(1 − ϕ12 L12 ). Porque as´ı se justificar´ıa el rezago 13 como el producto L × L12 = L13 . Entonces una primera decisi´on es una estructura SARMA (1, q)(1, q12)[12]. En los rezagos de media m´ovil se observan rezagos de media m´ovil 2,12,14. Con lo cual una opci´on podr´ıa ser colocar (1, 2)(1, 1)[12] porque as´ı se justificar´ıa el rezago 14 que proviene de 2+12 = 14. Es el modelo siguiente. (1 − ϕ1 L)(1 − ϕ12 L12 )Yt = (1 +

2 X

θj Lj )(1 + θ12 L12 )t ,

(6.39)

j=1

La estimaci´on de los modelos SARMA La estimaci´on del modelo (6.39), Yt ∼ ARM A(1, 2)(1, 1)[12], se hace con los comandos siguientes. mod4 = arima(r,order=c(1,0,2), seasonal = list(order = c(1, 0, 1), period = 12)) Este modelo tiene residuos dados por r4 = mod4$residuals. El an´alisis de incorrelaci´on de e´ stos se hace de la manera explicada en el curso, en los talleres y en las notas de clase. Igual con los pron´osticos.

Una aplicaci´on de los modelos ARMA estacionales, SARMA Estos modelos son una alternativa a los modelos con componentes estacionales y de tendencia. Por ejemplo, un modelo con tendencia lineal y posible componente estacional, St , seg´un el modelo de descomposici´on est´a dado por Yt = β0 + β1 t + St + t , t ∼ ARM A(p, q). La alternativa es un modelo de la forma Yt = β0 + β1 t + t ,

(6.40)

118 t ∼ SARM A(p, q)(ps, qs).

(6.41)

La idea es que en el modelo (6.41) el ARMA estacional modela la parte St + t del modelo (6.40).

CAP´ITULO

Ra´ıces Unitarias y Tendencias Estocasticas ´ (ARIMA)

7.1. Introducci´on En este cap´ıtulo se presenta en modelo alterno al de Descomposici´on con errores ARMA, conocido como modelo ARIMA-SARIMA o´ modelo de Box-Jenkins. El modelo ARIMA puede tomarse como an´alogo del modelo con tendencia, Yt = Tt + εt , y el modelo SARIMA como an´alogo del modelo con tendencia y estacionalidad, Yt = Tt + St + εt . Primero se introducen los modelos ARIMA, luego las pruebas de ra´ız unitaria de Dickey-Fuller para detectar estos modelos. M´as adelante se introducen los modelos SARIMA y las pruebas de ra´ız unitaria estacional para detectarlos.

7.2. MOdelos ARIMA Observaci´on: Si en Yt = a+bt+εt , εt ∼ RB(0, σ 2 ) se toma ∆Yt = (1−L)Yt = Yt −Yt−1 , se obtiene Yt − Yt−1 = a + bt + εt − (a + b(t − 1) + εt−1 ) = b + εt − εt−1 = b + εt + θ1 εt−1 , 119

θ1 = −1,

120 luego, llamando Wt = ∆Yt y ηt = εt − εt−1 , se obtiene Wt = b + ηt. N´otese que ηt es MA(1) y por tanto Wt es un M A(1) con media diferente de cero, estacionario en covarianza. Luego diferenciar una serie no estacionaria, con tendencia lineal, puede producir una serie estacionaria en covarianza. Definici´on 7.2.1 (Ra´ız Unitaria Autoregresiva). Suponga un modelo ARM A(p, q) ϕ(L)pYt = θ(L)q εt ,

εt ∼ RB(0, σ 2 ),

tal que una de las p ra´ıces de la ecuaci´on ϕp (z) = 1 − ϕ1 z − · · · − ϕp z p = 0 es 1, entonces se dice que Yt tiene una ra´ız unitaria autorregresiva. En este caso el polinomio ϕp (L) factoriza como ϕp (L) = (1 − L)ϕp−1 (L) donde ϕp−1 es un polinomio de grado p − 1 y ϕp−1 (1 − L)Yt = θ(L)εt , luego si Yt ∼ ARM A(p, q) con una ra´ız unitaria entonces ∆Yt ∼ ARM A(p − 1, q). Ejemplo 7.2.1. Si ϕ(L) = (1 − 0.2L)(1 − L),

θ(L) = 1 − 0.3L

(1 − 0.2L)(1 − L)Yt = (1 − 0.2L)εt

(1 − 1.2L + 0.2L2 )Yt = εt − 0.3εt−1

Yt − 1.2Yt−1 + 0.2Yt−2 = εt − 0.3εt−1 . Entonces Yt ∼ ARM A(2, 1), y no es estacionario en covarianza pues 1 − 1.2z + 0.2z 2 = 0 tiene ra´ıces z1 = 1, z2 = 5 y s´olo una cumple |z| > 1. Pero Wt = ∆Yt = Yt − Yt−1 cumple (1 − 0.2L)Wt = (1 − 0.3)εt ∼ ARM A(1, 1) es estacionario en covarianza e invertible. Definici´on 7.2.2. La marcha aleatoria sin tendencia se define como un proceso Yt ∼ AR(1), de media cero, con coeficiente ϕ = 1, es decir Yt = Yt−1 + εt , t ∈ Z,

εt ∼ RB(0, σ 2 ).

(7.1)

La marcha aleatoria con tendencia se define como un proceso AR(1) de media diferente de cero, con coeficiente ϕ = 1, es decir Yt = µ + Yt−1 + εt , µ ∈ R. La marcha aleatoria sin tendencia cumple, para t ≥ 0, Yt = Y0 + Y0 y V ar(Yt ) = tσ 2 .

(7.2) Pt

i=0 εi . Luego E(Yt |Y0 )

121 Para la marcha aleatoria con tendencia Yt = tδ + Y0 +

t X

εi ,

(7.3)

i=0

E(Yt|Y0 ) = tδ + Y0 ,

(7.4)

V ar(Yt) = tσ 2 .

(7.5)

0 Y.sin

−8

−4

40 20

Y.con

Un ejemplo de trayectorias de marchas aletorias con y sin tendencia se muestran a continuaci´on en la Figura 7.1.

60 Time

(a) con tendencia

100

Time

(b) sin tendencia

Figura 7.1: Trayectorias de Marchas Aleatorias Obs´ervese que la marcha aleatoria con tendencia Yt , tiene una caracter´ıstica similar a una serie con tendencia lineal de la forma Xt = a + bt + εt , εt ∼ RB(0, σ 2), porque en promedio crecen en cada per´ıodo [t − 1, t], una cantidad constante, E(Yt − Yt−1 ) = µ y E(Xt − Xt−1 ) = b. Suponga Yt = µ + ϕYt−1 + εt , t = 1, 2, . . ., T , un proceso AR(1) con media diferente de cero. Los pron´osticos a j pasos a partir de T son YˆT (j) = µ + ϕj YT para j = 1, 2, . . ., h. En el caso ϕ = 1, se obtiene una marcha aleatoria con tendencia. En este caso YˆT (j) = µ + YT , j ≥ 1, lo cual se interpreta como que una marcha aleatoria con tendencia no puede pronosticarse. Su mejor pron´osticos es el u´ ltimo valor conocido. P Finalmente, en una marcha aleatoria con tendencia, se tiene Yt = Y0 + µt + tj=1 εj , y, por tanto, para k = 1, 2, . . ., Cov(Yt , Yt+k ) = E((Yt − E(Yt ))(Yt+k − E(Yt+k )))

122

= E(

t X j=1

= σ t,

εj

t+k X

εj ) =

j=1

t X t X

E(εj εi )

j=1 i=1

Como V ar(Yt ) = σ 2 t entonces p Corr(Yt , Yt+k ) = σ 2 t/ σ 2 t(σ 2 (t + k)). p = 1/ 1 + k/t

Para t grande tal que k/t sea muy peque˜no se puede aproximar y obtener que las autocorrelaciones en marchas aleatorias con tendencia son casi uno, ρ(k) ≈ 1. Como la tendencia no determina el valor de la autocorrelaci´on, la conclusi´on tambi´en es v´alida en el caso de marchas aleatorias sin tendencia. Las marchas aleatorias son series con fuertes autocorrelaciones. Definici´on 7.2.3. Si la serie Yt no estacionaria en covarianza es tal que ∆Yt es estacionario, se dice que Yt es integrada de orden 1, I(1). Si ∆2Yt = ∆(∆Yt ) = (1−L)2 Yt es estacionario en covarianza se dice que integrada de orden 2, I(2). La definici´on de proceso integrado de orden d ≥ 1 es inmediata. Definici´on 7.2.4. Una serie Yt , t ∈ Z se dice que tiene tendencia estoc´astica si es integrada de orden d = 1, 2, I(d). La marcha aleatoria con tendencia (7.2) es un ejemplo de tendencia estoc´astica, en contraposici´on a una seria con tendencia lineal Tt = a + bt, la cual se dice que tiene tendencia determin´ıstica. Una serie con tendencia estoc´astica presenta per´ıodos de duraci´on aleatoria durante los cuales crece linealmente. Es el caso de los precios de las acciones. En per´ıodos de tiempo el precio puede aumentar por efecto de la demanda. Pero llega un momento en el cual la demanda se satura y los precios empiezan nuevamente a oscilar, o a tener otro comportamiento. La implementaci´on de los modelos ARIMA empieza con series con tendencia estoc´astica. La metodolog´ıa consiste en diferenciar la series d veces hasta obtener otra serie estacionaria en covarianza, identificable por un modelo ARMA, y estimar los par´ametros de este modelo. Definici´on 7.2.5 (Proceso ARIMA(p,d,q)). Una serie de tiempo Yt , t ∈ Z sigue un modelo ARIM A(p, 1, q), con tendencia, si (1 − L)Yt sigue un proceso ARMA(p,q) estacionario en covarianza, con media diferente de cero. Es decir, si ϕ(L)(1 − L)Yt = δ + θ(L)εt , y ϕ(z) = 0 tiene las ra´ıces fuera del c´ırculo unitario.

εt ∼ RB(0, σ 2),

(7.6)

123 Una serie de tiempo sigue un modelo ARIM A(p, d, q) si ∆d Yt = (1 − L)d Yt ,

d = 1, 2, . . .

es un proceso ARMA(p,q) estacionario en covarianza. En la pr´actica los casos m´as importantes son d = 0, 1, ver Diebold [1999, p´ag. 210]. De acuerdo con la definici´on (7.2.3) de serie integrada si Yt ∼ ARIM A(p, 1, q) entonces es integrada. Ejemplo 7.2.2. Si Yt ∼ ARM A(1, 1) con E(Yt) = 0 entonces (1 − ϕL)Yt = (1 + θL)εt y se puede comprobar que las autocorrelaciones de Yt satisfacen: ρ(1) =

(1 + ϕθ)(ϕ + θ) 1 + θ2 + 2ϕθ

ρ(k) = ρ(k − 1),

k ≥ 2.

Si ϕ → 1 obtenemos un ARIM A(0, 1, 1) = IM A(1, 1) y ρ(1) →

(1 + θ)(1 + θ) = 1. 1 + θ2 + 2θ

Por lo tanto ρ(k) ≡ 1, ∀k ≥ 1. N´otese que el IMA(1,1) es similar a la marcha aleatoria sin tendencia. Esto d´a lugar a afirmar que “los procesos ARIM A(p, 1, q) se comportan como marchas aleatorias en ciertos aspectos”. Ejemplo 7.2.3. En Diebold [1999, p´ag. 12] se considere la serie del Producto Interno Bruto en EUA entre 1869-1933, con el per´ıodo 1934-1993 para comparaci´on. Se examinan dos posibles modelos para esta serie. 1. Modelo con tendencia determin´ıstica lineal y errores AR(2) Yt = β0 + β1 t + εt εt = ϕ1 εt−1 + ϕ2 εt−2 + ηt,

ηt ∼ RB(0, σ 2 ).

2. Modelo con tendencia aleatoria Yt ∼ ARIM A(1, 1, 0) con tendencia (1 − ϕL)(1 − L)Yt = δ + εt

Zt = δ + ϕZt−1 + εt ,

Zt = (1 − L)Yt

La primera diferencia de Yt sigue un modelo AR(1) con media diferente de cero.

124 Entre 1932 y 1933 hubo una grave recesi´on en EUA y los pron´osticos con base en 1933 se hacen en una posici´on por debajo de la tendencia

Figura 7.2: Pron´osticos de la Serie PNB-USA.

El modelo 1 produce mejores pron´osticos que siguen la tendencia, en cambio, el modelo 2 no, pues subestima los valores. Ejemplo 7.2.4. An´alisis de la Tasa de cambio: USD/YEN, el precio de 1 US D´olar en Yen, ver Diebold [1999, p´ag. 221]. La serie es mensual desde 1973-01 hasta 1996-07. Por ejemplo, Yt = 107 significa que en el mes t se pagaban 107 Yen por 1 USD. El an´alisis siguiente es con el logaritmo de esta tasa. En lo que sigue y denota el logaritmo de la tasa. (b)

1985

1990

1995

0.0

(c)

(d)

ACF

0.00

1.0 Lag

−0.05

1980

0.5

Time

−0.10

1975

0.0 0.2 0.4 0.6 0.8 1.0

ACF 1980

0.05

1975

1985

1990

1995

Time

1.5

2.0

1.5

2.0

−0.2 0.0 0.2 0.4 0.6 0.8 1.0

4.4 4.6 4.8 5.0 5.2 5.4 5.6

(a)

0.0

0.5

1.0 Lag

0.2 0.1 −0.1

0.0

Partial ACF

0.3

(e)

0.5

1.0

1.5

2.0

Lag

Figura 7.3: (a) : Serie log(Precio) del Yen en USD, (b) Su fac, (c): La serie diferenciada, (d),(e): fac y facp de la serie diferenciada

125 La Figura 7.3 (b) sugiere un proceso integrado y las Figuras (d) y (e) sugieren un ARMA para la serie diferenciada. Utilizando la funci´on z = autoarmafit(dy), de la librer´ıa timsac, con la serie diferenciada, dy = diff(y,1,1), se obtiene un modelo ARMA(3,2). Utilizando la funci´on mod = arima(y,order=c(3,1,2)) se obtienen los resultados de la Tabla 7.1. Finalmente, utilizando la prueba Ljung-Box, se comprueba Tabla 7.1: Par´ametros del modelo ARIMA(3,1,2) parametros sd.dev est t ar1 0.18 0.06 2.91 ar2 -0.93 0.01 -77.71 ar3 0.37 0.06 6.17 ma1 0.20 0.02 11.32 ma2 0.99 0.03 32.43

que los residuos de este modelo son Ruido Blanco. El resutado de la prueba es X-squared = X-squared = 15.7889, df = 26, p-value = 0.9411,

5.0 4.5

log(tasa)

5.5

y no rechaza la hip´otesis de ruido blanco ya que el valor p es mayor de 0.05. Este modelo es v´alido para pronosticar el logaritmo de la tasa de cambio, a corto plazo. Hay que anotar que la funci´on arima() de R tiene varias versiones, y que no es lo mismo estimar el modelo ARIMA(3,1,2) con la variable log(Yt), que estimar el modelo ARMA(3,2) con la variable ∆ log(Yt ). Los pron´osticos a 19 per´ıodos se muestran en la Figura 7.4, comparados con los pron´osticos de un modelo con tendencia lineal y error estructural AR(2).

160

180

200

220

240

260

280

mes

Figura 7.4: Pron´osticos usd/yen con ARIMA(3,1,2)(continua) y Tendencia Lineal+AR(2) (punteada) Ejemplo 7.2.5. Considerando Yt la serie del logaritmo de la tasa de cambio USD/Libra,

126

−2

diff(Y, 1, 1)

−14 −10 −6

−2

es decir, el log del precio de 1 libra esterlina en usd, mensual, en el per´ıodo 01/80-12/88, con 478 observaciones. La trayectoria parece una marcha aleatoria sin tendencia, como se observa en la Figura 7.6.

100

Time

(a) Serie USD/Libra

100

Time

(b) Primera Diferencia de la Serie USD/Libra

Figura 7.5: Serie USD por Libra La primera diferencia parece estacionaria, casi ruido blanco. La fac de la primera diferencia confirma que es ruido blanco, obteniendo un modelo de marcha aleatoria con tendencia, dado por (1 − L)Yt = µ + εt ,

µ = −0.0009.

Este modelo se denomina Modelo de Black-Scholes o´ Modelo Log-Normal para el precio de un activo. Ejemplo 7.2.6. Sea Yt la serie del empleo en Canad´a de periodicidad trimestral en el per´ıodo 01/1962-04/1993, desestacionalizada. Se propuso un modelo AR(2). Pero usando la funci´on auto.arima, se obtiene Yt ∼ ARIM A(1, 1, 0) sin intercepto, es decir, (1 − ϕL)(1 − L)Yt = εt ,

εt ∼ RB(0, σ 2),

(7.7)

con ϕˆ = 0.4598, s.e. = 0.0758, σ ˆ = 2.068, BIC = 491.22, AIC = 485.41, M AP E = 1.0453. Si se ajusta un modelo AR(2) a la serie, es decir, Yt = µ + ϕ1 Yt−1 + ϕ2 Yt−2 + εt ,

(7.8)

se obtiene ϕ ˆ1 = 1.4505, ϕˆ2 = −0.4763, σ ˆ = 2.022, µ ˆ = 97.498, AIC = 493.57, BIC = 505.22, M AP E = 1.0648.

127 La conclusi´on es que ser´ıa preferible utilizar el modelo (7.7) por tener un mejor ajuste. Pero es posible utilizar otra herramienta para discriminar entre los posibles modelos. La prueba de hip´otesis Diebold-Mariano, en la funci´on dm.test de la librer´ıa forecast compara las medias de los errores de pron´ostico a h pasos, dentro de la muestra, para dos modelos. La hip´otesis nula es que las medias son iguales y la alterna que no lo son. Es posible utilizar hip´otesis de una cola, con la especificaci´on de que el primer modelo tiene menor error medio que el segundo. Se obtiene el valor-p para la prueba de dos colas: p-value = 0.8412. Luego no rechaza la hip´otesis nula y ambos modelos tienen la misma capacidad predictiva. #### programa en R library(forecast) # Para la funci´ on arima y la prueba dm.test # Canadian employment index, seasonally adjusted, 1961:1-1994:4 # 136 obs E = read.table("CAEMP.DAT", header = TRUE) y = ts(E$caemp, frequency = 4, start = c(1961,01), end = c(1994,04)) f1 = arima(y, order=c(2,0,0)) f2 = arima(y, order=c(1,1,0)) accuracy(f1) accuracy(f2) dm.test(residuals(f1), residuals(f2), h=1) ----------resultado data: residuals(f1) residuals(f2) DM = -0.2004, Forecast horizon = 1, Loss function power = 2, p-value = 0.8412 alternative hypothesis: two.sided

7.3. Ra´ıces Unitarias, Estimaci´on y Prueba de Hip´otesis Minimos Cuadrados Ordinarios (MCO) y Ra´ıces Unitarias La discusi´on en general es para un Yt ∼ ARIM A(p, 1, q) pero se puede utilizar una marcha aleatoria sin tendencia. Suponga un proceso AR(1), Yt = ϕYt−1 + εt , en el cual ϕ = 1 pero ϕ se estima usando MCO. El estimador de MCO de ϕ, ϕ, ˆ tiene dos propiedades 1. Superconsistencia

128 a) En el caso ϕ = 1 el estimador de MCO ϕˆT tiene la propiedad de que d

T (ϕˆT − 1) −→ Z,

T → ∞,

(7.9) d

donde Z es una variable alatoria no degenerada, y el s´ımbolo −→ denota la convergencia en distribuci´on. b) En el caso |ϕ| < 1, el estimador de MCO, ϕˆT cumple √

T (ϕˆT − ϕ) −→ Z

T → ∞,

(7.10)

donde Z es variable aleatoria no degenerada. √ como T < T , se dice que la convergencia en el caso (7.9) es m´as r´apida que en el caso (7.10), y a esto se denomina superconsistencia, es decir, el estimador de MCO, ϕˆT de una ra´ız unitaria es superconsistente. 2. Sesgo de ϕˆ Si ϕˆT es el estimador de MCO de ϕ entonces E(ϕˆT ) < ϕ y el sesgo es mayor cuando ϕ = 1. El sesgo es ϕ − ϕˆT y crece si se considera tendencia. Aunque ϕˆT converge a ϕ, cuando T → ∞, el sesgo puede ser apreciable en muestras no muy grandes.

La fac y facp muestrales en presencia de Ra´ıces Unitarias En un proceso Yt con ra´ız unitaria, por ejemplo, ARIM A(p, 1, q) la fac muestral converge a cero con mucha lentitud y la facp muestral presenta un valor cercano a 1 en el rezago k = 1 y el resto son aproximadamente cero. Por ejemplo la Figura 7.6 siguiente de la fac y la facp de la tasa USD/Libra. Series Y

0.6 −0.2

0.2

Partial ACF

0.6 0.2 −0.2

ACF

1.0

Series Y

10 Lag

Lag

Figura 7.6: FAC y FACP Muestral de un ARIMA(p,1,q).

129 Hay que tener una precauci´on con relaci´on a tomar la decisi´on de diferenciar la serie para buscar un ARMA cuando se observan las caracter´ısticas anteriores. Y es que es necesario aplicar una prueba de hip´otesis para determinar si existe una ra´ız unitaria. 1. Diebold [1999, p´ag. 221]: “Si no hay ra´ız unitaria es conveniente utilizar modelos de niveles (componentes determin´ısticas) y es adecuando diferenciar s´olo en el caso de ra´ız unitaria; si la diferenciaci´on es inadecuada puede ser da˜nina, incluso asint´oticamente”. 2. Cuando una serie tiene ra´ız unitaria la serie es no estacionaria y los estimadores MCO no se distribuyen normal. 3. Soares and Medeiros [2008, p´ag. 4]: “La mayor´ıa de los art´ıculos sobre pron´osticos de demanda de energia toman diferencias sin hacer una prueba previa de ra´ız unitaria. Esto es un error grande cuando la tendencia es deterministica, tomar diferencias introduce una componente MA no invertible el cual causa problemas serios de estimaci´on.”

7.3.1.

Prueba Dickey-Fuller (DF)

La prueba Dickey-Fuller se basa en asumir que la serie se puede aproximar por un proceso AR(1) con tres variantes: media cero, media diferente de cero y tendencia lineal. Inicialmente se asume que Yt sigue un modelo AR(1) y se procede a transformar el modelo de la siguiente manera. Yt = ϕ1 Yt−1 + εt , Yt − Yt−1 = (ϕ1 − 1)Yt−1 + εt , ∆Yt = ρYt−1 + εt .

donde ρ = ϕ1 − 1. La existencia de una ra´ız unitaria equivale a ϕ1 = 1, es decir, a ρ = 0. 1. Prueba DF para el caso 1: suponiendo que Yt ∼ AR(1) con media cero, entonces ∆Yt = ρYt−1 + εt .

(7.11)

La hip´otesis nula es H0 : ρ = 0 versus la alterna Ha : ρ < 0. El estad´ıstico de la prueba se denota por τ y su distribuci´on bajo H0 permite calcular los valores cr´ıticos, de tal forma que el criterio de rechazo es τˆ < τ0.05 , con τˆ es valor calculado del estad´ıstico. En R la prueba DF se implementa mediante la librer´ıa urca por medio de la funci´on ur.df(y,type="none",lags=0).

130 2. Prueba DF para el caso 2: suponiendo que Yt ∼ AR(1) con media diferente de cero, entonces ∆Yt = α + ρYt−1 + εt , (7.12) con la misma hip´otesis. En R es ur.df(y,type="drift",lags=0). 3. Prueba DF para el caso 3: suponiendo que Yt ∼ AR(1) con tendencia lineal, entonces ∆Yt = α + β t + ρYt−1 + εt ,

(7.13)

con la misma hip´otesis. En R es ur.df(y,type="trend",lags=0).

7.3.2.

Prueba Dickey-Fuller Aumentada

La prueba aumentada de Dickey-Fuller no es solamente una prueba sino que requiere una estrategia de an´alisis para su aplicaci´on. Como se˜nalan Elder and Kennedy [2001, pag. 139]: “ Un ingrediente crucial en esta prueba, que no se reconoce bien en los libros de texto, es que se requiere una estrategia de prueba, en oposici´on al simple c´alculo del estad´ıstico de la prueba. Esta estrategia es necesaria para determinar si un intercepto, un intercepto m´as una tendencia con el tiempo, o ninguna de las dos anteriores deber´ıa inclu´ırse al correr la regresi´on para la prueba de raiz unitaria. Inclu´ır demasiados regresores puede generar una p´erdida de potencia, mientras que no inclu´ır suficientes puede generar resultados sesgados hacia el no rechazo de la hip´otesis nula...inclu´ır intercepto, o´ intercepto m´as tendencia, es necesario para permitir una representaci´on de la hip´otesis alterna que pueda competir contra la hip´otesis nula de ra´ız unitaria.” La estrategia para la prueba Dickey-Fuller Aumentada consiste en determinar cu´al de los siguientes tres casos determina una mejor aproximaci´on a la serie original, Yt . N´otese que los casos a examinar dependen del order autorregresivo p, por lo que la b´usqueda para por ejemplo, p = 1, 2, 3, requiere 3 × 3 casos. La decisi´on se toma con base en el menor AIC. Caso 1 Suponiendo que Yt ∼ AR(p) con media cero, (“none”) entonces Yt =

p X

ϕj Yt−j + εt = ϕ1 Yt−1 +

j=1

Defina ρ1 =

p X

ϕj Yt−j + εt .

j=2

p X j=1

ϕj , ρ i = −

p X j=i

ϕj , i = 2, . . . , p,

(7.14)

131 entonces con ρ = ρ1 − 1, la ecuaci´on (7.14) se transforma en ∆Yt = ρYt−1 +

p X

ρj ∆Yt−j+1 + εt .

(7.15)

j=2

Si hay una ra´ız unitaria se cumple 1 − ϕ1 − ϕ2 − · · · − ϕp = 0, es decir, ρ1 = 1. En este caso se tiene que el modelo (7.15) equivale a Yt − Yt−1 = ∆Yt =

p X j=2

p−1 X

ρj (Yt−j+1 − Yt−j ) + εt ρj+1 ∆Yt−j + εt

(7.16)

j=1

es decir un AR(p − 1) en la variable Zt = ∆Yt . La hip´otesis nula es H0 : ρ = 0 y la alterna Ha : ρ < 0. El estad´ıstico DF, τˆ tiene la misma distribuci´on asint´otica que el estad´ıstico DF para ρ = 0, ver caso 1 (7.11). ((As´ı, los resultados del proceso AR(1) se generalizan asint´oticamente en forma directa a procesos de orden superior)). Diebold [1999, p´ag. 128] Caso 2 Suponiendo Yt ∼ AR(p) con media diferente de cero (“drift”). Con la misma notaci´on del caso 1, Yt − µ =

p X j=1

ϕj (Yt−j − µ) + εt

Yt = α + ϕ1 Yt−1 + ∆Yt = α + ρYt−1 +

p X

j=2 p X

ρj (Yt−j+1 − Yt−j ) + εt

ρj ∆Yt−j+1 + εt .

(7.17)

j=2

donde α = µ(1 −

j=1

ϕj ), con los ρj definidos como en el caso anterior.

La hip´otesis nula es H0 : ρ = 0 y la alterna Ha : ρ < 0. Bajo H0 : ρ = 0 el t´ermino P α se anula porque pj=1 ϕj = 1. La distribuci´on asint´otica del estad´ıstico DF es igual a la del caso AR(1) con media, ver caso 2 (7.12). Caso 3 Modelo AR(p) con tendencia lineal (“trend”). En este caso se define Yt = a + bt +

p X j=1

ϕj Yt−j − a − b(t − j) + εt

(7.18)

132 que se puede reordenar de la siguiente forma ∆Yt = k1 + k2 t + ρYt−1 +

p X

ρj ∆Yt−j+1 + εt .

(7.19)

j=2

donde k1 = a 1 − k2 = b 1 − ρ1 =

p X i=1

p X

i=1 p X i=1

ϕi ϕi

p X

iϕi,

i=1

ϕi , ρ = ρ1 − 1.

La hip´otesis nula es H0 : ρ = 0 y la alterna Ha : ρ < 0. Bajo H0 : ρ = 0 se tiene que Pp k2 = 0, k1 = b i=1 iϕi y el estad´ıstico DF tiene la misma distribuci´on asint´otica del estad´ıstico en el caso AR(1), ver caso 3 (7.13). No rechazar la hip´otesis nula de ra´ız unitaria no necesariamente significa que se deba asumir que existe, debido a la baja potencia en estas pruebas. Como se˜nala Diebold [1999, p´ag 220]: “Las pruebas de ra´ız unitaria tienen problemas de potencia y tama˜no de muestra. De potencia porque las hip´otesis alternas son muy cercanas a la hip´otesis nula”. La presencia de cambios estructurales en la serie se sabe que disminuye la potencia de las pruebas ADF, pero no se desarrollar´a este tema. Ejemplo 7.3.1. Continuando con el Ejemplo 7.2.4, del logaritmo de la tasa de cambio USD/Yen, prcio de 1 usd en yens, en donde se encontr´o que la serie ln Yt puede ser modelada por un IMA(1), es decir por un ARIMA(0,1,1). A continuaci´on se realizan las prueba DF y DF aumentada. El resultado es que en ambos casos no se rechaza la hip´otesis de ra´ız unitaria, por lo que el modelo ARIMA(0,1,1) queda justificado. La Prueba Dickey-Fuller. Se realiza la prueba de ra´ız unitaria tomando en cuenta la tendencia. Es decir, se aplica el caso 3: “trend”. ## analisis con pruebas dickey-fuller de yen/usd library(forecast) library(urca) library(fUnitRoots) uno = read.table("yen_usd.dat",header=T,stringsAsFactors=F)

133 attach(uno) y = ts(y,frequency=12, start=c(1973,01), end = c(1996,07)) ly = log(y) ## prueba df caso 3 = trend, con libreria urca df.trend = ur.df(y = ly, lags = 0, type = ’trend’ ) summary(df.trend) --------------- Resultados (salida de R) Test regression trend Call: lm(formula = z.diff ˜ z.lag.1 + 1 + tt) Value of test-statistic is: -1.7764 2.568 1.5811 Critical values for test statistics: 1pct 5pct 10pct tau3 -3.98 -3.42 -3.13 phi2 6.15 4.71 4.05 phi3 8.34 6.30 5.36 La prueba DF para el caso 3 de un AR(1) con tendencia lineal, a un nivel de significaci´on de 5 %. El estad´ıstico DF observado fu´e -1.7764 mayor que el valor cr´ıtico -3.42, luego no rechaza la hip´otesis nula de raiz unitaria. Aunque hay que tener en cuenta la baja potencia de la prueba. En conclusi´on, es factible un modelo ARIMA(p,1,q) para el logaritmo de la tasa de cambio yen-usd. La Prueba Dickey-Fuller Aumentada. A continuaci´on se presenta una implementaci´on de la estrategia de an´alisis mencionada en Elder and Kennedy [2001, pag. 139]. El objetivo es determinar cu´al caso de los tres posibles con modelos AR(p) es el que mejor aproxima la serie para as´ı lograr una mayor potencia en la prueba. Para esto se utiliza la librer´ıa dynlm que extiende la funci´on de regresi´on lineal lm() para series de tiempo, permitiendo inclu´ır valores rezagados Yt−j , con el comando L(y,j), tendencia lineal a + bt con el comando trend(y). Por ejemplo, para programar el modelo de la ecuaci´on (7.19), en la pag. 132, ∆Yt = k1 + k2 t + ρYt−1 +

p X

ρj ∆Yt−j+1 + εt ,

(7.20)

j=2

con p = 3, queda ∆Yt = k1 + k2 t + ρYt−1 + ρ2 ∆Yt−1 + εt , se programa reg.trend1 = dynlm(dly ˜ trend(ly) + L(ly,1) + L(dly, 1)).

134 y la correspondiente prueba DF aumentada se programa con la funci´on ur.df() de la librer´ıa urca con: ur.df(y = ly, lags = 2, type = ’trend’). ## estrategia de regresiones para la df aumentada require(dynlm) reg.drift0 reg.drift1 reg.drift2 reg.drift3

= = = =

dynlm(dly dynlm(dly dynlm(dly dynlm(dly

˜ ˜ ˜ ˜

L(ly,1)) L(ly,1) + L(dly, 1)) L(ly,1) + L(dly, 1) + L(dly, 2)) L(ly,1) + L(dly, 1) + L(dly, 2) + L(dly, 3))

(c(AIC(reg.drift0),AIC(reg.drift1),AIC(reg.drift2),AIC(reg.drift3))) reg.trend0 = reg.trend1 = reg.trend2 = reg.trend3 = + L(dly, 3))

dynlm(dly dynlm(dly dynlm(dly dynlm(dly

˜ ˜ ˜ ˜

trend(ly) trend(ly) trend(ly) trend(ly)

+ + + +

L(ly,1)) L(ly,1) + L(dly, 1)) L(ly,1) + L(dly, 1) + L(dly, 2)) L(ly,1) + L(dly, 1) + L(dly, 2)

(c(AIC(reg.trend0),AIC(reg.trend1),AIC(reg.trend2),AIC(reg.trend3))) ## se detecta el caso reg.trend1 como el de menor aic df.trend = ur.df(y = ly, lags = 2, type = ’trend’ ) summary(df.trend) -----resultados de los valores AIC para los 8 modelos -1202.376 -1240.019 -1236.218 -1232.318 -1203.215 -1245.147 -1240.917 -1237.595 -----resultados de la prueba DF aumentada Test regression trend lm(formula = z.diff ˜ z.lag.1 + 1 + tt + z.diff.lag) Value of test-statistic is: -2.58 2.9804 3.448 Critical values for test statistics:

135 1pct 5pct 10pct tau3 -3.98 -3.42 -3.13 phi2 6.15 4.71 4.05 phi3 8.34 6.30 5.36 A partir de este u´ ltimo resultado se concluye que no se rechaza la hip´otesis de ra´ız unitaria ya que el valor observado de la prueba -2.58 es mayor que el valor cr´ıtico al nivel de 5 %, -3.42. Luego, es v´alido diferenciar la serie. En el Ejemplo 7.2.4 se ajust´o el modelo ARIMA(3,1,2) que produce buenos pronosticos del logaritmo de la tasa.

136

CAP´ITULO

Ra´ıces Unitarias Estacionales y Estacionalidad Estocastica ´ (SARIMA)

8.1. Modelos SARIMA Si la serie Yt tiene una componente estacional con per´ıodo s es posible eliminarla diferenciando una o´ dos veces, con un rezago de orden s, es decir, transformando Yt a Wt = (1 − Ls )D Yt = ∆D s Yt ,

D = 1, 2,

(8.1)

esperando que la serie Wt sea estacionaria en covarianza y as´ı proceder a buscar una estructura ARM A.

8.1.1.

Modelo SARIMA

Sin embargo, si hay estacionalidad, pueden existir estructuras ARIM A intra y entre los s per´ıodos. Las ARIMA intra se refieren a modelos de la forma ϕp(L)∆d Yt = θq (L)εt,

εt ∼ RB(0, σ 2 ),

(8.2)

con d = 0, 1, 2. Los ARIMA entre se refieren a modelos de la forma s Φps (Ls )∆D s Yt = Θqs (L )εt ,

137

εt ∼ RB(0, σ 2 ),

(8.3)

138 con D = 0, 1, 2. Y se define el modelo SARIM A(p, d, q)(ps, D, qs)s como s ϕp (L)Φps (Ls )∆d ∆D s Yt = θq (L)Θqs (L )εt ,

(8.4)

d s D La transformaci´on (filtro lineal) Xt = ∆d ∆D s Yt = (1 − L) (1 − L ) Yt es la que elimina la tendencia y la estacionalidad de Yt , dejando una estructura ARM A(p + sps , q + sqs ). En el caso d = D = 0 se obtiene un modelo SARM A(p, q)(ps, qs ), (6.36), pag. 115.

La idea es que la serie filtrada Xt en el modelo ϕp (L)ΦP (Ls )Xt = θq (L)ΘQ (Ls )εt,

εt ∼ RB(0, σ 2).

tiene un papel similar o equivalente al residuo estructural εt en el modelo Yt = a + b t +

s−1 X

δj Ij (t) + εt .

(8.5)

j=1

En (8.4) al diferenciar la serie Yt se eliminan la tendencia y estacionalidad aleatorias. Pero diferenciar en (8.5) para eliminarlas puede ser un procedimiento incorrecto ya que estas componentes son b´asicas para calcular los pron´osticos. Este hecho llama la atenci´on sobre la necesidad de realizar pruebas de hip´otesis que permitar decidir entre ambos modelos cuando hay evidencia de estacionalidad. Ejemplo 8.1.1. Suponga un modelo SARIM A(1, 1, 0)(1, 1, 0)12, es decir, p = 1, d = 1, q = 0, ps = 1, D = 1, qs = 0, s = 12, dado por (1 − ϕ1 L)(1 − ϕ12 L12 )(1 − L)(1 − L12 )Yt = εt , Colocando Xt = (1 − L)(1 − L12 )Yt = Yt − Yt−1 − Yt−12 + Yt−13 se obtiene, de manera equivalente, (1 − ϕ1 L)(1 − ϕ12 L12 )Xt = εt . Cuando se desarrollan los polinomios autoregresivos del modelo se obtiene (1 − ϕ1 L − ϕ12 L12 + ϕ1 ϕ12 L13 )Xt = εt , o´ tambi´en, Xt = ϕ1 Xt−1 − ϕ12 Xt−12 + ϕ1 ϕ12 Xt−13 + εt . Note que este modelo se puede considerar como un AR(13) con restricciones sobre ciertos coeficientes, que deben asumirse cero. El par´ametro ϕ1 ϕ12 se estima como un solo valor ϕ13 y puede ser no significativo. En el caso que sea no significativo el modelo se dice sin interacci´on, y se expresa como (1 − ϕ1 L − ϕ12 L12 )Xt = εt .

139 Ejemplo 8.1.2. Consideremos el modelo del ejemplo anterior, con par´ametros ϕ1 = 0.8, ϕ12 = −0.137, σ 2 = 2. El siguiente c´odigo en R utiliza la funci´on sarima.Sim() de la librer´ıa CombMSC para simular trayectorias de modelos SARIMA y as´ı visualizar mejor esta clase de modelos. library(CombMSC) library(forecast) ti=0.8; ti4=-0.137; sigma= sqrt(2) Y = sarima.Sim(n = 22, period =12, model=list(order=c(1,1,0),ar = ti, ma = NULL, sd = sigma), seasonal=list(order=c(1,1,0),ar=ti4,ma = NULL))

−10

−5

ts.plot(Y) auto.arima(Y)

Time

Figura 8.1: Trayectoria de un modelo SARIMA(1,1,0)(1,1,0)[12] Ejemplo 8.1.3. Modelo “Air-passengers” Se define como un modelo multiplicativo de la forma SARIM A(0, 1, 1)(0, 1, 1)12 dado por (1 − L)(1 − L12 )Yt = (1 + θ1 )(1 + θ12 L12 )εt ,

|θ1 | < 1, |θ12 | < 1.

(8.6)

140 Es un modelo utilizado con frecuencia para modelar series con tendencia lineal y componente estacional. Por ejemplo la serie The classic Box Jenkins airline data. Monthly totals of international airline passengers, 1949 to 1960. . Nota sobre la identificaci´on de este modelo Al utilizar Wt = (1 − L)(1 − L12 )Yt y examinar la fac de Wt se observa que Wt = (1 + θ1 L)(1 + θ12 L12 )εt es un M A(13). Se puede comprobar que la fac de Wt tiene valores diferentes de cero s´olo en los rezagos 1,11,12,13. La funci´on de autocorrelaci´on, fac, de Wt es

2 ρ0 = (1 + θ12 )(1 + θ12 ) 2 ρ1 = θ1 (1 + θ12 )

ρ11 = θ1 θ12 ρ12 = θ12 (1 + θ12 )

Suponga el modelo air-passengers con par´ametros θ1 = −0.377, θ12 = −0.572, σ 2 = 0.0014, el siguiente c´odigo simula una trayectoria de este modelo.

#---------------------------------------------------------------# simular ti=-0.377; ti4=-0.572; sigma= sqrt(0.0014) Y = sarima.Sim(n = 12, period =12, model=list(order=c(0,1,1),ma = ti, ar = NULL, sd = sigma), seasonal=list(order=c(0,1,1),ma=ti4,ar = NULL)) # graficar ts.plot(Y) # identificar modelo = auto.arima(Y) summary(modelo) # estimar arima(Y,order=c(0,1,1),seasonal=list(order=c(0,1,1),period=12))

141

Time

Figura 8.2: Trayectoria de un modelo SARIMA(0,1,1)(0,1,1)[12] Ejemplo 8.1.4. Sea un proceso SARIM A(1, 1, 1)(1, 1, 1)4 definido por (1 − ϕ1 L)(1 − ϕ4 L4 )(1 − L)(1 − L4 )Yt = (1 − θ1 L)(1 − θ4 L4 )εt se simula con el c´odigo siguiente. Es interesante examinar si en R es equivalente estimar el modelo SARIMA con la serie Yt , y el modelo ARMA con la serie diferenciada, Wt = (1 − L)(1 − L4 )Yt . library(CombMSC) library(forecast) fi=-0.8; fi4=-0.37; ti=-0.64; ti4=-0.513; sigma= sqrt(0.014) Y = sarima.Sim(n = 60, period =4, model=list(order=c(1,0,1), ar = fi, ma = ti, sd = sigma), seasonal=list(order=c(1,1,1),ar=fi4,ma = ti4)) ts.plot(Y) auto.arima(Y) arima(Y,order=c(1,1,1),seasonal=list(order=c(1,1,1),period=4))

142 W = diff(diff(Y,1,1),4,1) acf(W) pacf(W) arima(W,order=c(1,0,1),seasonal=list(order=c(1,0,1),period=4))

8.2.

Pruebas de Ra´ız Unitaria Estacional

A partir de la definici´on de los modelos SARIMA hay dos posibles modelos para una serie con componente estacional: 1. Modelos con componente determin´ıstica, es decir, el modelo de descomposici´on. 2. Modelos no estacionarios, integrados estacionalmente, es decir, SARIMA(p,d,q)(P,D,Q)[s], con d 6= 0 y D 6= 0. Sin embargo, un modelo de la forma SARIMA(p,0,q)(P,0,Q)[s], equivalente a un ARMA(p+P,q+Q), podr´ıa ser tanto estacionario como mostrar caracter´ısticas estacionales. Por tanto, un tercer modelo es el anterior, que se denomina un modelo “estacionario estacional”. El modelo con componente determin´ıstica no es estacionario, pero si se elimina la componente de tendencia, el proceso filtrado resulta la suma de un proceso determin´ıstico peri´odico y otro estacionario en covarianza. Concretamente, Yt − Tˆt = Sˆt + εˆt . Este modelo tambi´en se denominar´a “estacionario estacional”. Se necesitan t´ecnicas para discriminar cu´al modelo aplicar. Aunque en el curso se han implementado criterios tales como la validaci´on cruzada y la comparaci´on del MAPE para escoger el mejor modelo para pronosticar, en la teor´ıa de series de tiempo se ha dado mucho desarrollo en las pruebas de ra´ıces unitarias, por ejemplo, las pruebas para detectar ra´ıces unitarias estacionales, como las pruebas de Hylleberg y otros, indicada por HEGY, y la prueba Canova-Hansen, indicada por CH, que se exponen a continuaci´on, con el prop´osito de utilizarlas como herramientas para escoger el modelo m´as adecuado. Para introducir la idea de las pruebas de ra´ız unitaria estacional hay que recordar que una serie Yt ∼ ARM A(p, q) se dice integrada, o´ con una ra´ız unitaria si se cumple que el polinomio autorregresivo ϕp(z) tiene una ra´ız igual a 1. Es decir, si se cumple que ϕp(1) = 0. En este caso el polinomio se puede factorizar como ϕp (z) = (1 − z)ϕp−1 (z). Definici´on 8.2.1. Una serie de tiempo Yt ∼ ARM A(p, q) se dice integrada estacional, de per´ıodo s con s par, si se cumple que algunas o´ todas las ra´ıces de z s = 1 son ra´ıces del polinomio autorregresivo ϕp (z).

143 Las ra´ıces de z s = 1, para por ejemplo s = 4, 12, se denominan las ra´ıces de la unidad, y son s n´umeros complejos. Para el caso s = 4 se tiene que las cuatro ra´ıces de la unidad √ son 1, −1, i, −i, donde i = −1 es la unidad imaginaria. Para s = 12 las ra´ıces son √ √ √ √ ±1, ±i,−(1 ± i 3)/2, (1 ± i 3)/2, −( 3 ± i)/2, ( 3 ± i)/2. Con la representaci´on z = reiθ donde r = |z| y θ es el a´ ngulo que forma |z| con el eje X, se pueden identificar las cuatro ra´ıces unitarias con los a´ ngulos θ = 2πj/4, j = 0, 1, 2, 3, ya que r = 1 en todos los casos, pues las ra´ıces unitarias tienen m´odulo 1. Estos a´ ngulos se denominan las “frecuencias estacionales”.

8.2.1.

La prueba HEGY

La prueba HEGY es una generalizaci´on de la prueba aumentada de Dickey-Fuller. Permite decidir cu´ales de las frecuencias estacionales corresponden a ra´ıces estacionales significativas. En caso de detectarse se concluir´a que la serie en cuesti´on est´a integrada estacionalmente integrada. Esta secci´on describe la idea b´asica de la prueba, que utiliza una regresi´on lineal m´ultiple. Los estad´ısticos para detectar las ra´ıces unitarias estacionales son tipo t-Student y F, correspondientes a los estad´ısticos para pruebas de hip´otesis sobre los par´ametros de la regresi´on. Lo que sigue est´a tomado del art´ıculo de Hylleberg et al. [1990]. Inicialmente se asume que la serie Yt est´a generada por un proceso AR(p), es decir, se cumple ϕ(L)Yt = εt , para εt ruido blanco. Se asume que est´an dados p n´umeros complejos θj ∈ C, j = 1, . . . , p, en los cuales ϕ(θj ) es un n´umero finito. Estos p n´umeros son las p ra´ıces estacionales. El valor p es el per´ıodo, indicado anteriormente por s. Adicionalmente se definen las cantidades auxiliares: δk (z) = 1 − z/θk , p Y ∆(z) = δk (z), k=1

λk = ϕ(θk )/

δj (θk ).

(8.7)

j6=k

Entonces se cumple, por una identidad (debida a Lagrange), que ϕ(z) =

p X λk ∆(z)(1 + δk (z)) k=1

δk (z)

+ ∆(z)ϕ∗ (z),

(8.8)

144 donde ϕ∗ (z) es un polinomio posiblemente infinito o´ racional. N´otese que, por la definici´on (8.7) se cumple ϕ(θk ) = 0 ⇔ λk = 0

(8.9)

por lo tanto, para determinar si θk es ra´ız estacional se debe comprobar λk = 0. Pero estas λk se pueden asimilar a los coeficientes en una regresi´on lineal m´ultiple que se define a partir de la identidad (8.8). Y por tanto, probar λk = 0 se puede hacer, en principio, mediante una prueba t-Student de significaci´on de par´ametros, o´ con una prueba F. Por ejemplo, para el caso ϕ(z) = 1 − z 4 , donde θ1 = 1, θ2 = −1,θ3 = i,θ1 = −1, la identidad (8.8) se transforma en ϕ(z) = λ1 z(1 + z)(1 + z 2 ) + λ2 (−z)(1 − z)(1 + z 2 )

+λ3 (−iz)(1 − z)(1 + z)(1 − iz) + λ4 (iz)(1 − z)(1 + z)(1 + iz) +ϕ∗ (z)(1 − z 4 ).

(8.10)

¯ 3 . Si se definen nuevas constantes Pero se tiene que λ3 , λ4 ∈ C, y adem´as, λ4 = λ π1 , π2 , π3, π4 mediante las relaciones: π1 = −λ1 , π2 = −λ2 , 2λ3 = −π3 + iπ4 , 2λ4 = −π3 − iπ4 entonces la identidad (8.10) queda ϕ(z) = −π1 z(1 + z + z 2 + z 3 ) − π2 (−z)(1 − z + z 2 − z 3 ) −(π4 + π3 z)(−z)(1 − z 2 ) + ϕ∗ (z)(1 − z 4 ).

(8.11)

Finalmente, a partir de ϕ(L)Yt = εt , se reemplaza ϕ(L) por la correspondiente expresi´on obtenida al reemplar la expresi´on para ϕ(z) en (8.11). Si se definen las variables X1,t = (1 + L + L2 + L3 )Yt , X2,t = −(1 − L + L2 − L3 )Yt ,

X3,t = −(1 − L2 )Yt ,

X4,t = (1 − L4 )Yt , se obtiene la ecuaci´on

ϕ∗ (L)X4,t = π1 X1,t + π2 X2,t + π3 X3,t + εt .

(8.12)

Colocando ϕ∗ (L) = 1, la ecuaci´on (8.12) es una regresi´on lineal m´ultiple. A partir de esta ecuaci´on se plantea la prueba de hip´otesis HEGY para el caso de per´ıodo s = 4.

145

Descripci´on de la Prueba La hip´otesis en la prueba HEGY, para el caso s = 4, se compone de tres hip´otesis, una para cada ra´ız unitaria estacional: −1, ±i. N´otese que la ra´ız z = 1 corresponde a una ra´ız unitaria para tendencia aleatoria. H0 : ϕ(1) = 0 ⇔ π1 = 0, Ha : ϕ(1) > 0 ⇔ π1 < 0, H0 : ϕ(−1) = 0 ⇔ π2 = 0, Ha : ϕ(2) > 0 ⇔ π2 < 0, H0 : |ϕ(i)| = 0 ⇔ π3 = π4 = 0, Ha : no(H0 ).

En el caso s = 12, hay 7 hip´otesis, dos para z = 1, −1, y una para cada pareja de ra´ıces unitarias estacionales conjugadas. Los estad´ısticos de las pruebas son t-Student para z = 1, −1, y F para ±i. Las decisiones se toman con base en los valores p correspondientes. N´otese que no rechazar una de las hip´otesis nulas equivale a aceptar que es una ra´ız unitaria estacional. En este caso hay que tener en cuenta la potencia de la prueba. Recordando que la potencia se refiere a la probabilidad: Prob(rechazar H0 | Ha es cierta), una baja potencia significa que la prueba no es capaz de detectar la alterna cuando e´ sta es cierta, o tambi´en, que la probabilidad de no rechazar H0 cuando es falsa, es alta. En caso de no rechazar una de las hip´otesis nula, cabr´ıa esperar que el modelo adecuado para la series sea de tipo SARIMA integrado, es decir, despu´es de diferenciar la serie se obtiene un proceso ARMA estacionario. En caso de rechazar todas las hip´otesis nulas la serie no tiene ra´ıces unitarias estacionales y entonces se clasifica como una serie “estacionaria estacional”, es decir, se puede modelar como un cierto ARMA con o´ rdenes altos, o´ tambi´en como un modelo de descomposici´on, como se vi´o al comienzo del curso.

Implementaci´on en R La prueba HEGY est´a implementada en la funci´on HEGY.test, de la librer´ıa uroot, (1 ). Al colocar res = HEGY.test(y,itsd,regvar,selectlags), la serie debe estar en la variable y, y debe haberse declarado como objeto “ts”. El par´ametros itsd se refiere a si se desea inclu´ır en el modelo un intercepto, tendencia lineal y componente estacional con 1

La librer´ıa uroot est´a en http://r-forge.r-project.org/bin/windows/contrib/latest/, en un archivo zip con el cual se instala directamente

146 variables indicadoras. El modelo se define en el caso afirmativo como (1 − L4 )(Yt − (β0 + β1 t +

3 X

δj Ij (t))) = εt .

(8.13)

j=1

Colocando itsd=c(1,1,c(1,2)) se especifica este modelo. En caso de que no se desee inclu´ır estas variables se coloca itsd=c(0,0,c(0)). El par´ametro regvar se coloca regvar=0 si no se incluyen variables ex´ogenas en el modelo. En caso contrario se incluye una matriz con los valores de las variables. Por ejemplo, se podr´ıa inclu´ır el primer rezago de y, Yt−1 . El par´ametro seleclags se coloca por defecto con los valores selectlags = list(mode = "bic", Pmax = 12). Ejemplo 8.2.1. Retomando el Ejemplo (5.1.1) en la pag. 80, para modelaci´on de la serie de producci´on de cemento Portland, trimestral se ajust´o un modelo de la forma Yt = β0 + β1 t +

3 X

δj It (t) + εt ,

(8.14)

j=1

εt ∼ AR(7).

(8.15)

500

−200

1000

s1 + e1

1500

200

400

2000

En la Figura 8.3 se puede ver la serie con y sin tendencia. Se observa que la componente estacional no es constante, por lo que es posible que se tenga un modelo SARIMA en lugar del modelo de componentes (8.14).

1960

1970

1980 Time

(a) serie con tendencia

1990

1960

1970

1980

1990

Time

(b) serie sin tendencia

Figura 8.3: Serie de Producci´on Trimestral de Cemento, Australia Con el fin establecer cu´al modelo es m´as adecuado se procede a realizar la prueba HEGY. Los comandos en R son

147 # prueba hegy res = HEGY.test(wts=y, itsd=c(0,0,c(0)), regvar=0, selectlags=list(mode="signf", Pmax=NULL)) res -------------resultados Null hypothesis: Unit root. Alternative hypothesis: Stationarity. ------------------------------------HEGY statistics: Stat. p-value tpi_1 3.245 0.1 tpi_2 -1.020 0.1 Fpi_3:4 2.254 0.1 Fpi_2:4 1.854 NA Fpi_1:4 4.101 NA

El resultado anterior muestra que no se rechaza la hip´otesis nula de ra´ıces unitarias en 1, −1, ±i, indicadas por tpi_1, tpi_2, Fpi_3:4. La u´ ltima corresponde a las ra´ıces en las frecuencias π/2, 3π/2, la primera corresponde a la frecuencia cero, y es una ra´ız unitaria para la tendencia, similar a la prueba Dickey-Fuller. Luego, es apropiado ajustar un modelo SARIMA a la serie. Al examinar la fac de la serie diferenciada wt = (1 − L)(1 − L4 )Yt , para Yt la serie excluyendo los u´ ltimos 8 datos para realizar validaci´on cruzada, se observa posibles modelos ARMA, por ejemplo, con rezagos 1, 12 y 13 significativos, seg´un la Figura 8.4.

100

0.0 −0.4

−0.2

ACF

200

0.2

148

Partial ACF

−0.3

−200

−0.1

−100

0.1 0.2 0.3

yi.d

Lag

1960

1965

1970

1975

1980

1985

1990

Time

Lag

(a) wt = (1 − L)(1 − L4 )Yt

(b) fac y facp de wt

Figura 8.4: Serie Diferenciada de Producci´on Trimestral de Cemento

Mediante ensayo y error, empezando con un modelo SARIMA(1,1,1)(1,1,1)[4], hasta conseguir residuos ruido blanco se llega al modelo SARIMA(3,1,2)(1,1,2)[4].

Tabla 8.1: Par´ametros estimados de SARIMA(3,1,2)(1,1,2)[4]

ar1 ar2 ar3 ma1 ma2 sar1 sma1 sma2

parametros 0.55 0.57 -0.31 -0.73 -0.23 0.90 -1.56 0.67

sd.dev 0.30 0.20 0.12 0.32 0.31 0.25 0.32 0.18

est t 1.85 2.82 -2.60 -2.28 -0.73 3.56 -4.83 3.80

La prueba Ljung-Box para los residuos de este modelo arroja el resultado X-squared = 34.672, df = 26, p-value = 0.1189, por lo que puede aceptarse este modelo para pronosticar. Los pron´osticos son el modelo SARIMA versus los obtenidos con el modelo componentes (8.14) se pueden observar en la Figura 8.5.

149

obs

1900

estr+ar(7)

1300

1400

1500

1600

1700

1800

sarima

Figura 8.5: Pron´osticos a 8 trimestres de la Producci´on Trimestral de Cemento

Como conclusi´on, se llega a que el SARIMA se detecta como modelo correcto por la prueba HEGY. Pero el MAPE de sus pron´osticos, a 8 trimestres, es 7.89, mayor que el del modelo de Descomposici´on con errores AR(7), que es 4.64. Sin embargo, el modelo SARIMA tiene mejor ajuste dentro de la muestra. En la Figura 8.6 se puede comparar tanto los ajustes dentro de la muestra como los pron´osticos. Cu´al modelo escoger?. Desde el punto de vista de la validaci´on cruzada a 8 trimestres, el modelo de Descomposici´on. Desde el punto de vista del ajuste dentro de la muestra, el modelo SARIMA.

500

1000

1500

2000

150

100

150

100

150

500

1000

1500

2000

c(t, tt)

50 c(t, tt)

Figura 8.6: Comparaci´on de los Ajustes y Pron´osticos. Panel superior: Modelo Descomposicion + AR, Panel inferior: SARIMA

8.2.2.

La prueba Canova-Hansen

Los modelos de Descomponsici´on con errores ARMA y los SARIMA se consideran, inicialmente, equivalentes. Pero, como se˜nalan Canova and Hansen [1995, pag. 238],

“Es dif´ıcil saber apriori cu´al posibilidad [modelo] produce la mejor descripci´on de la serie. Algunas series muestran cambios en los patrones estacionales, por ejemplo, el consumo de energ´ıa el´ectrica, las series de consumo e ingresos, la serie del producto interno bruto”.

La prueba CH se basa en esta idea. La inestabilidad estructural de la componente estacional se toma como una caracter´ıstica equivalente a la existencia de una ra´ız unitaria estacional y opuesta a estacionario estacional.

151

Descripci´on de la Prueba El modelo que asume la prueba es un modelo sin tendencia y con componente estacional, descrita mediante funciones trigonom´etricas, con el primer rezago de y como variable ex´ogena. Concretamente, la ecuaci´on (29) de la p´agina 243 en el art´ıculo de Canova and Hansen [1995], es Yt = bYt−1 +

s/2 X j=1

γj cos

2πjt s

γj∗ sen

2πjt , s

(8.16)

donde se asume s es par, t = 1, 2, . . .. N´otese que cuando j = s/2 entonces sin 2πjt = s sin(πt) ≡ 0. Luego, el n´umero de par´ametros en el modelo (8.16) es s − 1, el mismo del modelo con variables indicadoras. Los coeficientes se colocan en vectores γ = (γ1 , . . . , γs/2)0 , ∗ ∗ )0 . γ ∗ = (γ1∗ , . . . , γs/2−1 )0 γ ∗ = (γ1∗ , . . ., γs/2−1 La hip´otesis nula en la prueba CH es que los coeficientes γ, γ ∗ no cambian con t versus que var´ıan con t seg´un un modelo de marcha aleatoria, γ t = γ t−1 + ut , γ ∗t = γ ∗t−1 + vt , donde ut , vt son vectores iid, independientes de media cero. Esta hip´otesis alterna es la forma en la que se establece la presencia de ra´ıces unitarias estacionales, equivalente a inestabilidad estructural en la componente estacional. Las hip´otesis de la prueba se pueden escribir como sigue. H0 : Ha :

γ t ≡ γ, γ ∗t ≡ γ ∗ ,

γ t = γ t−1 + ut , γ ∗t = γ ∗t−1 + vt .

En Canova and Hansen [1995, pag. 240, sec. 2.2], se re-escribe la prueba con base en un par´ametro τ 2 , tal que H0 : τ 2 = 0, versus Ha : τ 2 > 0. El estad´ıstico de la prueba L, requiere resultados avanzados para su definici´on, por lo que se remite al art´ıculo original. Adem´as en el art´ıculo se proveen los valores cr´ıticos para varios niveles de significaci´on.

Implementaci´on en R La prueba Canova-Hansen est´a implementada en en la funci´on CH.test, que tiene los par´ametros siguientes. Al colocar res = CH.test(y,frec,f0,DetTr), la serie debe estar en la variable y, y debe haberse declarado como objeto “ts”. El par´ametro frec se refiere cu´ales frecuencias se prueban para detectar ra´ıces unitarias. N´otese que si s = 4 hay 2 frecuencias: −1, ±i, por lo que se coloca frec=c(1,1). En el caso s = 12 hay 6 frecuencias y se coloca frec=c(1,1,1,1,1,1). El par´ametro f 0 = 0, 1 se coloca igual a 1 si la variable Yt−1 se incluye como variable explicativa en el modelo, y se coloca igual a

152 0 si no se incluye. El par´ametro DetTr es una variable l´ogica con valores TRUE, FALSE. La opci´on TRUE indica que se incluye en el modelo una tendencia lineal. N´otese, sin embargo, que para aplicar la prueba Canova and Hansen [1995] es necesario eliminar la tendencia lineal de la serie Yt , en caso de existir. Para aplicar la prueba se prefiri´o filtrar previamente esta tendencia, colocando Yt − Tbt = Sbt + εbt , donde Tbt se estim´o mediante el filtro stl(). En consecuencia colocamos DetTr=FALSE. Ejemplo 8.2.2. Continuando el Ejemplo 8.2.1, para decidir entre los modelos de Descomposici´on con errores AR(7) y el SARIMA 1) Yt = β0 + β1 t +

3 X j=1

δj It (t) + εt .εt ∼ AR(7), 2) Yt = SARIM A(3, 1, 2)(1, 1, 2)[4].

La prueba se implementa con los comandos m1 = stl(y, s.window = ’per’, t.window = 50, t.jump = 1) s1 = m1$time.series[,1] t1 = m1$time.series[,2] e1 = m1$time.series[,3] y1 = s1+e1 ch.out1 = CH.test(wts=y1, frec=c(1,1,1), f0=0, DetTr=FALSE) ch.out1 --------------------Y los resultados son Canova & Hansen test ------ - ------ ---Null hypothesis: Stationarity. Alternative hypothesis: Unit root. Frequency of the tested cycles: pi/2 , pi , NA , L-statistic: 2.274 Critical values: 0.10 0.05 0.025 0.01 1.28 1.47 1.63 1.88

153 Como el estad´ıstico observado es 2.274, mayor que el valor cr´ıtico al nivel de 5 %, se rechaza la hip´otesis nula y se detecta una ra´ız unitaria estacional, que corresponde a inestabilidad estructural en la componente estacional. Por tanto, el modelo SARIMA es preferible para modelar la serie.

154

´ APENDICE

Datos de Series de Tiempo

A.1. Series con Tendencia y Estacionalidad A.1.1.

Ejemplo 1. Cinco Series

Para cargar los datos en R simplemente copiar cada bloque “structure” en el programa. Las series tiene componente estacional y/o tendencia.

"X1"

[PDF] Notas de Clase Series de Tiempo con R - Free Download PDF (2024)

FAQs

¿Qué son las series de tiempo según autores? ›

References