NOTAS DE CLASE DE SERIES DE TIEMPO Paul Castillo Bardalez y Fernando Perez 1 Universidad Nacional de Ingeniería Ciclo II 2005
1 Versión
preliminar por favor no circular, sugerencias son bienvenidas.
Índice general I Series de tiempo estacionarias
1
1. Introducción
1
1.1. Ecuaciones en diferencia . . . . . . . . . . . . . . . . . . . . . 1.2. Oper peradores de Rezago . . . . . . . . . . . . . . . . . . . . . .
2. Modelos ARIMA 2.1. 2.2. 2.2. 2.3. 2.4. 2.5. 2.5.
7
Introducción . . . . . . . . . . . . . . . . AR( p) p) . . . . . . . . . . . . . . Modelo Modeloss AR( Mode odelos MA(q) . . . . . . . . . . . . . . Mode odelos ARMA(p,q) . . . . . . . . . . . El teor teorem emaa de repr repres esen enta taci ción ón de Wold old .
. . . . . . . .. . .. .. ... . .. .. .. . . . . . .
. . . . .. .. . .
. . . . .
. . . 7 . .. 8 . . . 10 . . . 12 . . . 13
3. Predicción en modelos ARMA 3.1. 3.2. 3.2. 3.3. 3.3. 3.4. 3.4. 3.5. 3.5. 3.6. 3.6.
1 4
Conceptos básicos . . . . . . . . . . . Pred Predic icci ción ón con con inte interv rval alos os de con… con…an anza za Pred Predic icci ción ón con con pará paráme metr tros os esti estima mado doss . Calc Calcul ulan ando do el mejo mejorr pred predic icto torr line lineal al . Crit Criter erio ioss para para eval evalua uarr pred predic icci cion ones es . . El esta estadí díst stic icoo de Dieb Diebol oldd-Ma Mari rian anoo . .
4. Identi…caciòn de los modelos ARIMA
15 . . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
.. . . . . .. .. . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
.. . . . . .. .. . .
15 17 17 18 20 22
25
4.1. La metodo todollogía gía de Box Box y Jenk Jenkin inss . . . . . . . . . . . . . . . . 25
5. Estimación
29
5.1. La func unción ión de vero erosim similit litud . . . . . . . . . . . . . . . . . . . 30 5.2. Algoritmos Algoritmos numéricos numéricos de optimización optimización 1 . . . . . . . . . . . . . 33 5.2. 5.2.1. 1. Méto Método doss libr libres es de deri derivvadas adas . . . . . . . . . . . . . . . 34 1
Este capítulo se basa en el Libro .A pplied Computational Economics and Finance"de Mario J. Miranda y Paul L. Fackler v
Índice general I Series de tiempo estacionarias
1
1. Introducción
1
1.1. Ecuaciones en diferencia . . . . . . . . . . . . . . . . . . . . . 1.2. Oper peradores de Rezago . . . . . . . . . . . . . . . . . . . . . .
2. Modelos ARIMA 2.1. 2.2. 2.2. 2.3. 2.4. 2.5. 2.5.
7
Introducción . . . . . . . . . . . . . . . . AR( p) p) . . . . . . . . . . . . . . Modelo Modeloss AR( Mode odelos MA(q) . . . . . . . . . . . . . . Mode odelos ARMA(p,q) . . . . . . . . . . . El teor teorem emaa de repr repres esen enta taci ción ón de Wold old .
. . . . . . . .. . .. .. ... . .. .. .. . . . . . .
. . . . .. .. . .
. . . . .
. . . 7 . .. 8 . . . 10 . . . 12 . . . 13
3. Predicción en modelos ARMA 3.1. 3.2. 3.2. 3.3. 3.3. 3.4. 3.4. 3.5. 3.5. 3.6. 3.6.
1 4
Conceptos básicos . . . . . . . . . . . Pred Predic icci ción ón con con inte interv rval alos os de con… con…an anza za Pred Predic icci ción ón con con pará paráme metr tros os esti estima mado doss . Calc Calcul ulan ando do el mejo mejorr pred predic icto torr line lineal al . Crit Criter erio ioss para para eval evalua uarr pred predic icci cion ones es . . El esta estadí díst stic icoo de Dieb Diebol oldd-Ma Mari rian anoo . .
4. Identi…caciòn de los modelos ARIMA
15 . . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
.. . . . . .. .. . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
.. . . . . .. .. . .
15 17 17 18 20 22
25
4.1. La metodo todollogía gía de Box Box y Jenk Jenkin inss . . . . . . . . . . . . . . . . 25
5. Estimación
29
5.1. La func unción ión de vero erosim similit litud . . . . . . . . . . . . . . . . . . . 30 5.2. Algoritmos Algoritmos numéricos numéricos de optimización optimización 1 . . . . . . . . . . . . . 33 5.2. 5.2.1. 1. Méto Método doss libr libres es de deri derivvadas adas . . . . . . . . . . . . . . . 34 1
Este capítulo se basa en el Libro .A pplied Computational Economics and Finance"de Mario J. Miranda y Paul L. Fackler v
ÍNDICE GENERAL
vi
5.2. 5.2.2. 2. Méto Método doss basa basado doss en deri derivvadas adas . . . . . . . . . . . . . . 34
6. Aplicaciones de Laboratorio
37
6.1. 6.1. Iden Identi ti…c …cac ació iónn de los los proc proces esos os ARIM ARIMA A .. .. .. .. .. .. . 6.1.1. AR(1) . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.1.2. MA(1) . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.1.3. AR(2) . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.1.4. MA(2) . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.1.5 1.5. Comando ndos usu usuales ales en Evie Eview ws . . . . . . . . . . . . . . 6.1.6. 6.1.6. Apéndice: Apéndice: Simu Simulaci lación ón de procesos procesos ARMA ARMA en MATLA MATLAB B 6.1.7. AR(2) . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.1.8. MA(1) . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.1.9. Asignación . . . . . . . . . . . . . . . . . . . . . . . . 6.2. 6.2. Pred Predic icci ción ón y Exper Experim imen ento toss Comp Comput utac acio iona nale less . . . . . . . . . 6.2.1. Introdu oducción . . . . . . . . . . . . . . . . . . . . . . . 6.2.2 2.2. Linea ineam mient entos básico sicoss . . . . . . . . . . . . . . . . . . . 6.2.3. 6.2.3. Predicció Predicciónn y bandas bandas de con…anza con…anza en en un ARMA ( p; q ) 6.2.4. 6.2.4. Encontr Encontrando ando el mejor mejor predictor predictor lineal lineal en un AR( AR ( p) p) . . 6.2.5 2.5. Pred Prediicció cciónn en Evi Eviews ews . . . . . . . . . . . . . . . . . . . 6.2.6 2.6. Anál nálisis isis de indi indica cado dore ress . . . . . . . . . . . . . . . . . . 6.2.7. Ejercicios . . . . . . . . . . . . . . . . . . . . . . . . . 6.3. 6.3. Esti Estima maci ción ón por por Máxi Máxima ma Veros erosim imil ilit itud ud . . . . . . . . . . . . . 6.3.1. Introdu oducción . . . . . . . . . . . . . . . . . . . . . . . 6.3.2. 6.3.2. Algoritm Algoritmos os de Optimiza Optimización ción y Máxima Máxima Verosimil erosimilitud itud . 6.3.3 6.3.3.. Estim Estimaci ación ón por por Máxim Máximaa Vero Verosim simil ilitu itudd en Eview Eviewss . . . 6.3. 6.3.4. 4. Anál Anális isis is de la e…ci e…cien enci ciaa del del esti estima mado dorr . . . . . . . . . 6.3.5. Test de signi…cancia . . . . . . . . . . . . . . . . . . . 6.3.6 6.3.6.. Estim Estimaci ación ón por por Máxim Máximaa Vero Verosim simil ilitu itudd en Matlab Matlab . . . 6.3.7. Asignación . . . . . . . . . . . . . . . . . . . . . . . . 6.4. Ejercicios . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7. Vectores Autorregresivos Autorregresivos2 7.1. .1. 7.2. 7.2. 7.3. 7.3. 7.4. .4. 7.5. .5. 2
Esta Estacciona ionarrieda iedadd y moment entos . . . . Caus Causal alid idad ad en el sen sentido tido de Gran Grange gerr La repr repres esen enta taci ción ón en medi medias as móvi móvile less Funci nción impuls pulsoo resp respuuest esta . . . . . . Desc escompo ompossició iciónn de Varia rianza nza . . . . .
37 37 38 39 40 41 42 45 46 48 49 49 50 51 53 57 60 61 61 61 62 64 66 67 67 70 72
75 . . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
77 81 81 83 83
La mayoría del material discutido en estas notas es un resumen de los capítulos 10 y 11 de Hamilton (1994) y del libro de Amisano y Giannini "Topics in Structural VAR Econometrics"
ÍNDICE GENERAL
vii
7.6. Est Estimació aciónn de model odeloos VAR(p) . . . . . . . . . . . . . . . . . 85 7.7. 7.7. Iden Identi ti…c …cac ació iónn de mode modelo loss VAR(p AR(p)) . . . . . . . . . . . . . . . . 87
8. VAR Estructurales 8.1. Des Descom composi posici cióón de Chol holesk eski . . . . . 8.2. 8.2. Rest Restri ricc ccio ione ness con contemp tempor orán ánea eass . . . . 8.2.1. Bernanke y Mihov . . . . . . 8.2.2. Cushman y Zha . . . . . . . . 8.3. Varia riables bles no esta estaci cioonar narias ias . . . . . . 8.3 8.3.1. Rest estricc ricciiones nes de Larg argo Plaz Plazoo
89 . . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
91 92 92 95 96 97
9. Modelos VAR en la práctica
103
10. Modelos Var Estructurales en la Práctica
123
9.1. Introducción . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 9.2. Motivación . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 103 9.3. 9.3. El sistema sistema de Vecto Vectores res autorr autorregres egresivo ivos: s: Una Recap Recapitul itulació aciónn . 104 9.3.1. Estimación . . . . . . . . . . . . . . . . . . . . . . . . 104 104 9.3. 9.3.2. 2. Funci unción ón de impu impuls lsoo-re resp spue uest staa . . . . . . . . . . . . . . 106 106 9.3. 9.3.3. 3. Desc Descom ompo posi sici ción ón de Varia arianz nzaa . . . . . . . . . . . . . . . 106 106 3 9.4. 9.4. El Modelo de Stock Stock y Watson Watson (2001) (2001) . . . . . . . . . . . . . 10 1 07 9.5. Var en Eviews . . . . . . . . . . . . . . . . . . . . . . . . . . . 108 9.5.1. Estimación . . . . . . . . . . . . . . . . . . . . . . . . 108 108 9.5. 9.5.2. 2. Funci unción ón de impu impuls lsoo-re resp spue uest staa . . . . . . . . . . . . . . 111 111 9.5. 9.5.3. 3. Desc Descom ompo posi sici ción ón de Varia arianz nzaa . . . . . . . . . . . . . . . 112 112 9.5. 9.5.4. 4. Exte Extens nsió ión: n: Pred Predic icci ción ón en VAR . . . . . . . . . . . . . 113 113 4 9.6. VAR en MATLAB MATLAB . . . . . . . . . . . . . . . . . . . . . . . . 115 9.6.1. Estimación . . . . . . . . . . . . . . . . . . . . . . . . 115 115 9.6. 9.6.2. 2. Funci unción ón de impu impuls lsoo-re resp spue uest staa . . . . . . . . . . . . . . 120 120 10.1. Introducción . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123 10.2. El Modelo de Stock y Watson (2001) . . . . . . . . . . . . . . 123 10.3 10.3.. El modelo modelo con restr restricc iccion iones es de larg largoo plaz plazoo . . . . . . . . . . . 127 10.4. Midiendo los choques de Política Monetaria: Bernanke y Mihov (1998) (en progreso...) . . . . . . . . . . . . . . . . . . . . 131 10.4.1. Introdu oducción . . . . . . . . . . . . . . . . . . . . . . . 131 10.4.2. Identi…cación y estimación . . . . . . . . . . . . . . . . 133
3
Stock y Watson (2001), “Vector Autorregressions”, Journal of Economic Perspectives (JEP), Volumen 15, Número 4, Páginas 101-115. 4 La base de estos programas puede ser encontrada en "Econometric Toolbox for MATLAB"en http://www.spatial-econometrics.com
viii
ÍNDICE GENERAL
II Model odelos os de Seri Series es de Tiem iempo no Estac stacio iona nari rias as
137
11.Raíces Unitarias y Cointegración
141
10.5. Introducción . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139
11.1. Introducción . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141 11.2. Tests de Raíces Unitarias 5 . . . . . . . . . . . . . . . . . . . . 142 11.3. 11.3. El Test Test de Dic Dicke keyy – Full Fuller er aume aument ntado ado (ADF) (ADF) . . . . . . . . . 144 144 11.4. 11.4. El El Test Test de Dickey Dickey – Fulle Fullerr aumentado aumentado (ADF) (ADF) en Eviews Eviews . . . 146 11.5. 11.5. Raíces Raíces unitaria unitariass en presenci presenciaa de cambios cambios estru estructura cturales les . . . . 150
12.Cade 12.Cadenas nas de de Markov Markov y Modelo Modeloss de Cambio Cambioss de Régim Régimen en
153 153
12.1. Introducción y motivación . . . . . . . . . . . . . . . . . . . . 153 12.2. Cadenas de Markov . . . . . . . . . . . . . . . . . . . . . . . 154 12.2 12.2.1 .1.. De… De…ni nici cion ones es prin princi cipa pale less . . . . . . . . . . . . . . . . . 154 154 12.2.2. Distribuciones estacionarias . . . . . . . . . . . . . . . 156 12.2.3. Derivación de la probabilidades del estado estacionario 156 12.2 12.2.4 .4.. Est Estac acio iona nari ried edad ad Asin Asintó tóti tica ca . . . . . . . . . . . . . . . 157 157 12.2. 12.2.5. 5. Sim Simula ulació ciónn de de una una Caden Cadenaa de de Mark Markov ov . . . . . . . . . 157 157 12.3. Modelos de Markov Switching . . . . . . . . . . . . . . . . . . 159 12.3.1. Introdu oducción . . . . . . . . . . . . . . . . . . . . . . . 159 12.3.2. Estimación del modelo . . . . . . . . . . . . . . . . . . 160 12.3.3. Filtrado de probabilidades . . . . . . . . . . . . . . . . 162 12.3. 12.3.4. 4. Suav Suaviza izado do (Sm (Smoot oothin hing) g) de de proba probabi bilid lidade adess . . . . . . . 163 163 12.4. .4. Alguna gunass exten xtensi sion ones es . . . . . . . . . . . . . . . . . . . . . . . 165 12.4.1. Aproximación discreta . . . . . . . . . . . . . . . . . . 165 12.4 12.4.2 .2.. Mod Model elad adoo de pref prefer eren enci cias as . . . . . . . . . . . . . . . . 166 166 12.5. Referencias . . . . . . . . . . . . . . . . . . . . . . . . . . . . 167 12.6 12.6.. Apén Apéndi dice ce:: Prog Progra rama mass de de Mat Matla labb . . . . . . . . . . . . . . . . 168 168 12.6.1. ar1.m . . . . . . . . . . . . . . . . . . . . . . . . . . . 168 168 12.6. .6.2. Logli glike_a e_ar1.m 1.m . . . . . . . . . . . . . . . . . . . . . . 170 12.6.3. Filter.m . . . . . . . . . . . . . . . . . . . . . . . . . . 172 12.6.4. Smooth.m . . . . . . . . . . . . . . . . . . . . . . . . . 173 12.6.5. markov_ss.m . . . . . . . . . . . . . . . . . . . . . . . 174 12.6.6. simulate_markov.m 6 . . . . . . . . . . . . . . . . . . . 175 12.6.7. ejemplo_markov.m . . . . . . . . . . . . . . . . . . . . 175 12.6.8. tauchen.m tauchen.m 7 . . . . . . . . . . . . . . . . . . . . . . . . 176 12.6.9. consumo.m . . . . . . . . . . . . . . . . . . . . . . . . 177
5
Tomado de Chumacero (2001) Tomado del web site de Chris Edmond: http://pages.stern.nyu.edu/~cedmond/406/ 7 Tomado de www.eco.utexas.edu/~cooper/dynprog/tauchen.m 6
Parte I
Series de tiempo estacionarias
1
Capítulo 1
Introducción En este primer capítulo se presentan algunas herramientas básicas que son necesarias para el desarrollo futuro del curso. En primer lugar vamos a estudiar ecuaciones en diferencia lineales de orden p y operadores de rezago, que son de mucha utilidad para analizar modelos de series de tiempo estacionarias. El estudio de estas herramientas tiene como por objetivo familizar al estudiante con el manejo de la nomenclatura de los modelos de series de tiempo. El segundo grupo de herramientas busca proveer algunos algoritmos básicos para resolver problemas numéricos, en particular nos concentraremos en algoritmos para resolver problemas de optimización y para simular series de tiempo arti…ciales.
1.1. Ecuaciones en diferencia Las ecuaciones en diferencia relacionan variables discretas en el tiempo. Supongamos que y t y !t representan dos variables económicas, por ejemplo producto bruto interno y productividad. Una ecuación en diferencias de primer orden estará dada por: yt = y t1 + ! t
(1.1)
La ecuación anterior puede interpretarse de la siguiente manera: el producto bruto interno depende de su valor pasado y de la productivdad del periodo corriente. Nótese que la expressión (1.1) es una ecuación puesto que relaciona dos variables, una endógena y otra exógena. En este caso yt es la variable endógena, la variable para la cuál se busca una solución , mientras que ! t es la variable exógena. La solución de esta ecuación en diferencias de primer orden consiste en expresar yt como función únicamente de valores 1
2
CAPíTULO 1 INTRODUCCIÓN
presentes y pasados de !. La forma más sencilla de resover esta ecuación es sustituyendo recursivamente la misma expresión para periodos pasados, t 1, t 2, t 3, ... yt1 = yt2 + ! t1 yt2 = yt3 + ! t2
.. .
y1 = y 1 + ! 0
La solución, luego de sustituir recursivamente las ecuaciones anteriores, es la siguiente: yt = t+1 y1 + t ! 0 + t1 ! 1 + ::: + !t1 + ! t
(1.2)
Notese que la solución es un promedio ponderado de los valores presentes y pasados de la variable ! y de un valor inicial y1 : Cuando < 1 , los valores de ! más alejados en el pasado tienen menor importancia en la determinación de y t . Utilizando la ecuación (1.2) podemos de…nir multiplicadores dinámicos, estos muestran la respuesta de yt+ j a un cambio en la varaible ! t . Esta respuesta está dada por: @y t+ j = t @! t
(1.3)
Dependiendo del valor de , se pueden observar diferentes dinámicas en la respuesta de yt . Si 0 < < 1 el multiplicador dinámico decae geométricamente hacia cero. Mientras que si jj > 1 el sistema es explosivo. Cuando = 1, la variable endógena se convierte en la suma histórica de la variable exógena. yt = y1 + ! 0 + ! 1 + ::: + ! t1 + ! t
(1.4)
en este caso, @y@!t tj = 1, y se interpreta generalmente como que la serie tiene memoria in…nita puesto que los valores pasados de ! , independientemente de que tan distantes en el pasado se encuentren, mantienen su importancia en la determinación de y t .Este multiplicador dinámico es también conocido como función impulso respuesta. La ecuación en diferencias de orden 1 puede generalizarse al caso de una ecuación de orden p : +
yt = 1 yt1 + 2 yt2 + 3 yt3 ::: p yt p + ! t
(1.5)
3
1.1 ECUACIONES EN DIFERENCIA
La ecuación anterior puede resolverse de manera análoga al caso anterior, transformándola de forma adecuada utilizando notación matricial de la siguiente forma:
X t = F X t1 + t
Donde, las nuevas variables se de…nen com:
X t =
26 66 664
yt yt1 yt2
.. .
37 77 775
, F =
26 66 664
1 2 2 : : : p1 p 1
0
0 :::
0
0
0
1
0 :::
0
0
.. .
.. .
.. . :::
.. .
.. .
37 77 775
, t =
26 66 664
(1.6) !t 0 0
.. .
37 77 775
0 Esta notación equivale a expresar la ecuación en diferencias de orden p como un sistema de ecuaciones en diferencias de orde 1. Este sistema admite una solución similar al caso de una ecuación en diferencias de orden 1 de la siguiente forma: yt p
0
0
0 :::
1
0
X t = F t+1 X 1 + F t 0 + F 1t1 1 + ::: + F t1 + t
Una forma alternativa de resolver este sistema de ecuaciones en diferencia de orden 1 es transformándolo en un sistema recursivo de la forma: Z t = Z t1 + V t
donde se obtiene a partir de los valores propios de la matriz F, tal que: F = T T 1
y Z t = T 1 X t V t = T 1 t
Los valores propios de F pueden encontrarse a partir de las raíces del siguiente polinomio de grado p p + 1 p1 + 2 p2 ::: + p1 + p = 0
4
CAPíTULO 1 INTRODUCCIÓN
1.2. Operadores de Rezago En esta sección introducimos un nuevo operador matemático que se aplica a series de tiempo. Una serie de tiempo se de…ne como la colección de T un conjunto de observaciones indexadas por el tiempo, fyt gt=t t=t . Ejemplos de series de tiempo son: la tasa de interés para el periodo 1983-2005, el producto bruto interno : 1992-2005, etc. Ejemplos de series de tiempo sencillas son: Tendencia: y t = t Constante yt = c Ruido Blanco Gausiano: yt = t en donde t es un proceso estocástico, conjunto de variables aleatorias, independiente y distribuido normalmente con media cero y varianza constante: t N (0; 2 ). El operador de rezagos se aplica sobre series de tiempo, y transforma la serie de tiempo original en otra serie de tiempo. El operador se simboliza por la letra L, y tiene las siguientes propiedades: 0
yt = Lx t
! yt = xt1
yt = L k xt
! yt = xtk
yt = bL k xt
! yt = bxtk
yt = bL k xt
! yt = bxt+k
!
yt = b1 L1 + b2 L2 xt
yt = b 1 xt1 + b2 xt2
Una ecuación en diferencias de orden 1 puede escribirse de una manera conveniente utilizando el operador de rezagos de la siguiente manera: (1
L) yt = ! t
mientras que la ecuación en diferencias de orden p como:
1
1 L
2L2:: pL p
yt = ! t
Se puede simpli…car notación de…niendo el polinomio de rezago de orden p (L), como:
(L) = 1
1 L
2L2:: pL p
5
1.2 OPERADORES DE REZAGO
De esta forma, una ecuación en diferencias de orden p puede escribirse como: (L)yt = ! t
Para resolver ecuaciones en diferencia utilizando el operador de rezagos la siguiente identidad resulta de mucha utilidad: 1 = 1 + L + 2 L2 + 3 L3 :::::::::: 1 L
para todo jj < 1 . En el caso de la ecuación en diferencias de orden 1, la solución seria la siguiente: (1 yt =
L) yt = ! t
!t = 1 + L + 2 L2 + 3 L3 :::: ! t (1 L)
Con ello, la solución para la variable yt se puede escribir como sigue, después de aplicar las propiedades del operador de rezagos: k=
yt =
1
X k=0
k ! tk
Esta misma estratégia de solución puede generalizarse para el caso de una ecuación en diferencias de orden p. Utilizando la siguiente estratégia:
1
1 L
2L2:: pL p
yt = ! t
Notese, sin embargo que el polinomio de rezagos 1 1L 2 L2 :: p L p puede factorizarse convenientemente de la siguiente manera:
1
1 L
2L2:: pL p
= (1
1L) 1 2L2
Por ello, la solución puede escribirse como: yt =
(1
1L) (1
(1
1L) (1
1 2 L2 ) ::: (1
=
2L2)
!t 2 L2 ) ::: (1
::: (1
pL p)
pL p )
c1 c2 c p + +::: (1 1 L) (1 2 L2 ) (1 p L p )
6
CAPíTULO 1 INTRODUCCIÓN
Donde, cada uno de estos parámetros auxiliares pueden encontrarse de forma sencilla, utilizando el sistema de ecuaciones recursivas: ci =
(i
pi 1 1 ) (i 2 ) ::: (1
p)
(1.7)
Utilizando la notación previamente desarrollada la solución de una ecuación en diferencias de orden p estará dada por: yt =
c p c1 c2 ! t + ! t + ::: !t 2 (1 1 L) (1 2 L ) (1 p L p )
y de esta expresión se tiene que_
@y t+ j = c 1 j1 + c2 j2 + c3 j3 + :: + c p j p @! t
Capítulo 2
Modelos ARIMA 2.1. Introducción ´ En este ca´pitulo presentamos un conjunto de modelos basicos de series de tiempo estacionarias que servi´ran de base para el desarrollo de modelos mas ´ complejos en la segunda parte del curso. Una proceso estocastico, fyt g ´ se dice que es estacionario en el ´sentido debil si sus primeros momentos centrales son independientes del tiempo. E (yt ) =
E (yt E ((yt
)
2
= 2
) (yt j )) = j
Donde , 2 y j representan la media, varianza y covarianza de orden j del proceso esto´castico yt . Ejemplos de procesos esto´casticos simples son: un ruido blanco, " t , donde: E ("t ) = 0
E "2t = 2 E ("t "t+ j ) = 0
Dos funciones de mucha utilidad para caracterizar procesos esto ´casticos estacionarios son las funciones de autocorrelacíon parcial ( FAP) y simple 7
8
CAPíTULO 2 MODELOS ARIMA
(FAS). El FAS se de…ne como la correlacíon temporal entre yt y yt j , para todo j , y se de…ne: j =
j 0
De otro lado, el FAP denota la correlacíon entre y t y y t j una vez que se ha descontado el efecto intermedio de yt y yts para todo s < j . De…nimos la FAP como: yt = 11 yt1 + 22 yt1 + ::ss yts + t
donde ss , la correlacíon parcial de orden s . A continuacíon desarrollamos de manera detallada las caracte´risticas de los modelos autorregresivos, medias moviles ´ y ARMA models.
2.2. Modelos AR( p) Los modelos AR( p) tienen la siguiente representacíon: yt = 1 yt1 + 2 yt2 + 3 yt3 ::: p yt p + ! t
donde !t es un ruido blanco. El mismo modelo puede representarse utilizando el operador de rezago de la siguiente forma: (L)yt = ! t
donde el polinomino de rezagos (L) esta dado: (L) = (1
1L 2L2 3L3::: pL p)
Este polinomio de rezago puede factorizarse utilizando sus correspondientes ráices: (L) = (1
1L) (1 2L) :: (1 pL)
Utilizando la solucíon desarrollada en el ca´pitulo 1, podemos expresar el modelo AR( p) como un modelo MA (1). La forma MA del modelo AR(p) ofrece una serie de ventajas, en particular permite calcular todos los momentos de yt , los impulsos respuestas, y como veremos en el siguiente ca´pitulo, también obtener las correspondientes varianzas de los errores de prediccíon del modelo.
9
2.2 MODELOS AR(P ) k=
yt = +
1
X k=0
k ! tk
Donde: k = c 1 k1 +c2 k2 +c3 k3 +::+c p pk , y ci ha sido de…nido en ecuacíon (1.7 ). Utilizando, la presentacíon de medias ´moviles de…nida anteriormente permite calcular todos los segundos momentos de y t de forma sencilla: E (yt ) =
X X k=
E (yt
)
2
= 2!
1
2k
k=0
k=
E ((yt
) (yts
))
1
= 2!
k k+s
k=0
La condicíon necesaria y su…ciente que garantiza que los momentos de…nidos previamente sean …nitos es que las ráices del polinomio de rezago (L) se encuentren fuera del ´circulo unitario. ji j < 1:Un caso particular y ´facil de obtener es el caso de un modelo AR(1), en este caso, el modelo puede escribirse como: yt = 1 yt1 + ! t
En este caso, k = k1 . utilizando esta equivalencia, y las ´formulas derivadas para los segundos momentos de yt tenemos que: 0 =
2! 1 k1
s = s1 0
Una forma alternativa de encontrar las segundos momentos de un modelo AR(p) es utilizando las condiciones de Yule-Walker. Estas ecuaciones se obtienen multiplicando la ecuacíon AR(p) por yts , para todo s. Esto permite formar en un sistema ecuaciones lineales en s . Véamos un ejemplo con un modelo AR(2). yt = 1 yt1 + 2 yt2 + ! t
Multiplicando la ecuacíon anterior por yt ; aplicando el operador esperanza, bajo la condicíon de estacionariedad tenemos:
10
CAPíTULO 2 MODELOS ARIMA 0 = 1 1 + 2 2 + 2! 1 = 1 0 + 2 1 2 = 1 1 + 2 0 s = 1 s1 + 2 s2
El sistema anterior permite resolver para todos los valores de s en funcíon de los pa´rametros 1 ; 2 y : 2! : Una forma estratégia alternativa que se puede utilizar para invertir, obtener la representación de medias móviles, de un modelo autorregresivo es utilizar la siguiente identidad (1
1L 2L2 3L3::: pL p)(1 + 1L 2L2::::) = 1
por motivos de exposición asumamos que el orden del proceso autorregresivo es 2, es decir tenemos: (1
1L 2L2)(1 + 1L + 2L2::::) = 1
Multiplicando y agrupando ordenadamente este polinomio tenemos el siguiente sistema recursivo de ecuaciones: (1 + (1
1)L + (2 11 2) L2::: = 1
de donde se obtiene que: 1 = 1 2 = 1 1 + 2 , 3 = 2 1 + 1 2 ,
en general tenemos: s = s1 1 + 2 s2
2.3. Modelos MA(q) Los modelos de medias móviles tiene la siguiente representación general: yt = + ! t + 1 ! t1 + 2 ! t2 ::q !tq
11
2.3 MODELOS MA(Q)
donde, !t es un proceso ruido blanco. Utilizando operadores de rezago, un proceso MA(q) puede escribirse como: yt = + (L) ! t
donde, el polinomio de rezago, (L) se de…ne como: (L) = 1 + 1 L + 2 L2 ::q Lq
Note que todos los procesos MA(q) son estacionarios, independientemente del valor que tomen los parámetros i . Sin embargo, no todos los proceso MA(q) son invertibles, esto es, pueden transformarse en un proceso autorregresivo. De manera similar al caso de estacionariedad con modelos autorregresivos, en los modelos de medias móviles la condición de invertibilidad esta referida la valor de las raíces del polinomio de rezagos (L). En particular, el modelo presentado en la ecuación () será invertible si las q raíces del polinomio de rezago (L) se encuentran fuera del círculo unitario! t =
yt = (L) (1
1L) (1
yt 2 L) :: 1
qp L
donde j 1 j < 1: Los primeros y segundos de momentos pueden calcularse directamente aplicando el operador esperanza a yt , (yt )2 y (yt ) (yts ) de la siguiente manera. E (yt ) =
k=q
E (yt
X X )
2
= 2!
2k
k=0
k=q
E ((yt
) (yts
))
= 2!
k k+s
k=0
donde, E ((yt ) (yts )) = 0, para s > q:Para ganar más intuición sobre las propiedades de los modelos de médias móviles, consideremos el caso de un MA(1). Sus primer y segundos momentos estan dados por: E (yt ) =
E (yt
)
2
= 2! 1 + 21
12
CAPíTULO 2 MODELOS ARIMA E ((yt
) (yt1 )) = 2! (1 + 1)
Recuerde que en el caso de los modelos de médias móviles, la función de correlación parcial es cero para valores de s > q , y que la función de autocorrelación simple decae exponencialmente hacia cero.
2.4. Modelos ARMA(p,q) Un proceso ARMA(p,q) incluye tanto una parte autorregresiva como una de medias móviles. yt = c + 1 yt1 + 2 yt2 + 3 yt3 ::: p yt p + ! t + ! t + 1 ! t1 + 2 ! t2 ::q ! tq
Expresando el proceso en términos de operadores de rezagos tenemos que: (L)yt = c + (L) ! t
este proceso es estacionario si es que todas las raíces del polinomio de rezagos (L) se encuentran fuera del círculo unitario. En este caso, se puede encontrar la representación de medias móviles de este proceso, y expresar yt como: yt = + (L)! t
en donde, (L) = k=
tal que:
1
Xj k=0
(1
1 + 1 L + 2 L2 ::q Lq 1 L 2 L2 3 L3 ::: p L p )
j 1 y = 1 c :::
k <
1
2
3
p
Dado, los k , los momentos del proceso ARMA(p,q) pueden determinarse de la misma forma que para el caso de un AR(p). En caso que la serie que se este analizando no sea estacionario, posea algunas de las raíces del polinomio de rezagos de la parte autorregresiva dentro del círculo unitario, se puede modelar dicha serie utilizando un modelo ARMA(p,q) luego de transformar dicha serie en estacionaria. Para transformar una serie que posee d raíces unitarias en una seria estacionaria, hay que diferenciarla d veces. En donde, el operador diferencias se de…ne como:
2.5 EL TEOREMA DE REPRESENTACIÓN DE WOLD = (1
13
L)
Supongamos el siguiente proceso AR(3), tal que: (1
1L) (1 2L) (1 3L) yt
= ! t
en donde 2 = 3 = 1, entonces se puede decir que yt posee dos raíces unitarias.
1L) (1 L)2 yt = ! t De…niendo como z t = (1 L)2 yt = 2 yt , es claro que luego de tomar 2 (1
diferencias el proceso que queda, z t ; puede modelarse como un ARMA(1,0). Entonces se dice que yt sigue un proceso ARIMA (1,2,0). En general, un proceso ARIMA (p,d,q) es aquél que admite una representación ARMA(p,q) luego de tomarle d diferencias.
2.5. El teorema de representación de Wold Todo proceso que es estacionario en covarianzas puede escribirse como: yt = c +
1
X k=0
donde ,
Pj 1
k=0
k ! tk
j 1, y ! t es un ruido blanco.
k <
Proposition 1 i) (descomposición de Wold). Cualquier proceso estacionario en covarianzas, yt puede representarse de la forma: yt = c +
1
X k=0
donde: 0 = 1 y
1
P
2k <
k=0
k ! tk
1. El término ! t es un ruido blanco y repre-
senta el error de la predicción lineal de yt sobre sus valores rezagados ! t = y t E t (yt =yt1 ; yt2 :::yt p )
Capítulo 3
Predicción en modelos ARMA 3.1. Conceptos básicos Supongamos una serie de tiempo estacionaria y ergódica, fyt g con una representación de Wold dada por: yt = c +
1
X
k "tk
(3.1)
yt = c + "t + 1 "t1 + 2 "t2 ::::
Donde "t N (0; 2) , y además donde I t = fyt ; yt1 ::g denota el conjunto de información disponible en el periodo t, con primeros y segundos momentos …nitos, dados por : E (yt ) = c var(yt ) =
2
1
X
2k
k=0
La predicción de y ; h periodos en el futuro, y t+h ;condicional en la información I t se de…ne como yt+h=t , mientras que el error de predicción asociado como: t+h=t = y t+h
yt+h=t
El error cuadrático medio del error de predicción sirve de criterio de evaluación de la calidad de la predicción, y esta dado por: 15
16
CAPíTULO 3 PREDICCIÓN EN MODELOS ARMA
h
MSE (t+h=t ) = E 2t+h=t = E yt+h
i
yt+h=t 2
El mejor predictor de yt+h basado en la información I t , aquél que genera el error cuadrático medio mínimo, esta dado por la esperanza condicional de yt+h en el conjunto de información.
yt+h=t = E yt+h=I t
Prueba ( ver Hamilton páginas 72-73) Si f t g es un proceso ruido blanco distribuído independientemente, se tiene que E t+h=I t = 0 y que la predicción ópitma de y , E yt+h=I t es una función lineal de ft g
yt+h=t = c + h "t + h+1 "t1 ::
El error de predicción del mejor predictor lineal esta dado por: t+h=t = yt+h
yt+h=t
= "t+h + 1 "t+h1 ::: + h1 "t+1
(3.2) (3.3)
y su error cuadrático medio esta dado por:
MSE (t+h=t ) = 2 1 + 21 + :: + 2h1
Algunos comentarios:
E (t+h=t ) = 0 t+h=t no esta correlacionado con ningún elemento en I t
La forma de yt+h=t esta relacionada con la función de impulso respuesta MSE (t+h=t ) = var(t+h=t ) < var(yt ) lmh!1 yt+h=t = c lmh!1 MSE (t+h=t ) = var(yt )
3.2 PREDICCIÓN CON INTERVALOS DE CONFIANZA
17
3.2. Predicción con intervalos de con…anza si ft g se distribuye normalmente yt+h=I t
N (yt+h=t; 2
1 + 21 + :: + 2h1 )
En este caso un intervalo de con…anza de 95 por ciento para la predicción h periodos hacia adelante tiene la forma siguiente: :yt+h=I t
q
1;96
2 1 + 21 + :: + 2h1
3.3. Predicción con parámetros estimados En este caso el mejor predictor lineal con parámetros estimados esta dado por:
b b b b
= c + "t + "t1 :: yt+h=t h h+1
b
donde las variables con circun‡ejo denotan variables estimadas. El error de predicción con parámetros estimados esta dado por: t+h=t = yt+h
b
yt+h=t
b b
(3.4)
b
= "t+h + 1 "t+h1 ::: + h1 "t+1
b
(3.5)
Dado que j son variables aleatorias, el error cuadrático medio del error de predicción con parámetros estimados no coincide con el error de predicción del modelo teórico.
MSE (t+h=t ) = M SE (t+h=t ) = 2 1 + 21 + :: + 2h1
b
6
En este caso, el error cuadrático medio del error de predicción se puede descomponer entre el error de predicción teórico y otro termino que depende de la varianza del estimador de los parámetros.
h b b b i b b
MSE (t+h=t ) = E (yt+h
b
Operando adecuadamente,
+ y yt+h=t t+h=t
MSE (t+h=t ) = M SE (t+h=t ) + Y t E
b
yt+h=t )2
0
Y t0
18
CAPíTULO 3 PREDICCIÓN EN MODELOS ARMA
3.4. Calculando el mejor predictor lineal El mejor predictor lineal de yt+h=t puede calcularse de varias formas que son equivalentes. El algoritmo para calcular la predicción de un modelo AR(1) es bastante simple y metodologicamente se puede generalizar para calcular la predicción de un modelo ARMA, y para modelos multivariantes. Considere el siguiente modelo AR(1)
b
yt
c = (yt1 c) + t
donde " t N (0; 2 ) Para encontrar la predicción h periodos hacia adelante, se resuelve iterativamente hacia adelante: yt+h
c = h (yt c) + t+h + t+h1 + ::h1t+1
El mejor predictor lineal de yt+1 ,yt+2 ::yt+h se calculan usando la ley de proyecciones iteradas. yt+1=t
c = (yt c) yt+2=t c = yt+1=t c = 2 (yt c) yt+h=t c = (yt+h1 c) = h (yt c)
El correspondiente error de predicción esta dado por: t+1=t = y t+1 t+2=t = yt+2 t+h=t = y t+h
yt+1=t = t+1
yt+2=t = t+2 + t+1
yt+h=t = t+h + t+h1::::h1t+1
Las varianzas del error de predicción estas dadas por: var(t+1=t ) = 2
var(t+2=t ) = 2 1 + 2 var(t+h=t ) =
2
2
4
2(h 1)
1 + + ::
Tomando limites cuando h se aproxima a in…nito: lm yt+h=t = c
h
!1
=
1
2(h1) 1 2
3.4 CALCULANDO EL MEJOR PREDICTOR LINEAL lm var(t+h=t ) =
h
!1
El caso de AR(p)
(L) (yt
donde: (L) = 1
19
2 = var(yt ) 1 2
c) = t
1L 2L2:::: pL p
El algoritmo para la predicción de un modelo AR(p) es esencialmente el mismo que un AR(1), puesto que un modelo AR(p) se puede escribir como un modelo AR(1) utilizando la representación de estado espacio.
0B BB @
yt
1 C C CA c
yt1
c
:
yt p+1
c
=
0B B@
1 2 :: p 1 0 :: 0 0 :: 1 0
1C 0B CA BB@
yt1 yt2 : yt p
1 C C CA c c
+
c
o en términos de notación matricial, se tiene que:
0B BB @
t 0 : 0
1C CC A
t = F t1 + wt var(wt ) = w
Comenzando en el periodo t, y iterando hacia adelante h periodos, tenemos: t+h = F h t + wt+h + F wt+h1 :::::F h1 wt+1
Entonces el mejor predictor lineal basado en información hasta I t t+1=t = F t t+2=t = F 2 t t+h=t = F h t
El error de predicción asociado están dados por: wt+1=t = t+1
t+1=t = wt+1
20
CAPíTULO 3 PREDICCIÓN EN MODELOS ARMA wt+2=t = t+2
t+2=t = wt+2 + F wt+1
t+h=t = wt+h + F wt+h1 + :: + F hwt+1
wt+h=t = t+h
y el correspondiente matrices de varianza-covarianza del error de predicción. var(wt+1=t ) = w var(wt+2=t ) = w + F w F 0 h 1
var(wt+h=t ) =
X
F j w F j 0
j=0
Notese que: var(wt+h=t ) = w + Fvar(wt+h=t )F 0
3.5. Criterios para evaluar predicciones Para comparar la capacidad predictiva de un modelo con modelos alternativos se utilizan una serie de criterios que se basan en medidas de la varianza del error de predicción. Dado un horizonte de predicción, h; y n valores predichos, se tiene que: Raíz del error cuadrático medio
v ut X b T +n
RECM =
(yt
n
t=T +1
El error absoluto medio
T +n
MAE =
X bj (yt
t=T +1
yt )2
yt )
n
j
Media absoluta del porcentaje de error
T +n
MAPE = 100
X bj
t=T +1
Coe…ciente de desigualdad de Theil
(yt
yt ) =yt n
j
3.5 CRITERIOS PARA EVALUAR PREDICCIONES
s P b s P b s P T +n
t=T +1
UTHEIL =
T +n
t=T +1
21
(yt yt )2 n
(yt )2 n
T +n
t=T +1
(yt )2 n
Los dos primeros estadísticos del error de predicción dependen de la escala de la variable dependiente. Estos deberían usarse para comparar predicciones para la misma serie para diferentes modelos. Mas pequeño el error, mejor es la habilidad para predecir del modelo . Los dos últimos indicadores son independientes de la escala . La UTHEIL siempre entres cero y uno, cuando es cero indica un ajuste perfecto. El error cuadrático medio puede descomponerse como:
v ut X b X b b P b b T +h
t=T +1
yt )2
(yt
n
=
yt =n
2
y
+ sy
sy
2
+ 2(1
)s by sy
donde, yt =n, y ,sy , sy representan la media y la desviación estándar de yt y y t y es la correlación entre estas dos variables. . Se pueden de…nir las siguientes proporciones La proporción del sesgo:
b
Ps b P b s b P b (
T +n
(yt yt )2 n
t=T +1
La proporción de varianza
y)2
yt =n
2
sy
sy
T +n
(yt yt )2 n
t=T +1
La proporción de covarianza 2(1
)sy sy
s P b b T +n
t=T +1
(yt yt )2 n
22
CAPíTULO 3 PREDICCIÓN EN MODELOS ARMA
La proporción del sesgo nos dice cuan lejos de la media de la serie se encuentra la media de la proyección. La proporción de varianza, mide la diferencia entre la varianza de la proyección y la de la seria. Finalmente, la proporción de covarianzas mide los restantes componentes no sistemáticos del error de predicción. Si la proyección es buena, la proporción del sesgo y la varianza deberían ser pequeños.
3.6. El estadístico de Diebold-Mariano 1 Supongamos la serie fyt g para que se realizan dos predicciones yt+h=t y 2 yt+h=t obtenidas con dos modelos de predicción diferentes basados en información hasta el periodo t. Los errores de predicción de ambos modelos están dados por:
1t+1=t = y t+1
1 yt+1=t
2t+1=t = y t+1
2 yt+1=t
Para determinar que model predice mejor que otro, podemos realizar la siguiente prueba de hipótesis:
6
H o : E (L 1t+1=t ) = E (L 2t+2=t )
en contra de la hipótesis alternativa
H 1 : E (L 1t+1=t ) = E (L 2t+2=t )
donde,
L it+1=t = it+1=t
2
La prueba de Diebold-Mariano se basa en la diferencia de las funciones de perdida,
dt = L 1t+1=t
L 1t+2=t
Este estadístico se de…ne como: S =
donde:
b db p avar(d) = p d LRV=T
3.6 EL ESTADíSTICO DE DIEBOLD-MARIANO
b P P d =
1 T 0
23
dt
LRV = 0 + 2
1
j
j=1
Diebold y Mariano (1995) demuestran que bajo la hipótesis nula de igual capacidad predictiva, S
N (0; 1):
Capítulo 4
Identi…caciòn de los modelos ARIMA En las secciones anteriores se ha presentado en detalle las propiedades de los modelos ARMA (p,q), en particular hemos mostrado cuales son los primeros y segundos momentos de las series a partir de asumir que la serie bajo análisis, y t ; sigue un proceso en particular. En esta sección por el contrario, la idea es identi…car que modelo es el más apropiado para capturar la dinámica de la data, dados los primeros y segundos momentos muestrales de la serie. En particular, vamos a utilizar de manera recurrente las funciones de autocorrelación simple,FAS, y parcial, FAP para determinar los valores de p y q . Estas dos funciones estan relacionadas con los parámetros de los modelos ARMA(p,q), y contienen información relevante para el proceso de identi…cación.
4.1. La metodología de Box y Jenkins Esta metodología permite identi…car el proceso más apropopiado para la serie de tiempo yt ; utilizando un proceso iterativo de estimación y chequeo de las propiedades del modelo. Este proceso puede resumirse en los siguientes cuatro pasos: 1. Transformar la data hasta que el proceso se convierta en uno estacionario. Si una tendencia es evidente en los datos hay que removerla. Si la tendencia es determinística, esta se remueve estimando para la serie yt un polinomio determinìstico de la forma: yt = 0 + 1 t + 2 t2 ::s ts + "t
25
26CAPíTULO 4 IDENTIFICACIÒN DE LOS MODELOS ARIMA el residuo de esta regresión se considera como el proceso estacionario que se esta buscando identi…car. Por otro lado, si la tendencia no es determinística se procede a diferenciar la serie tantas veces sea necesario para convertirla en estacionaria. (L)d yt = (L)
Para determinar si se ha alcanzado estacionariedad se debe examinar las FAS y FAP de la serie transformada. Si esta es estacionaria, ambas debería aproximarse a cero rápidamente. 2. Hacer un supuesto inicial sobre los ordenes p y q del modelo ARMA(p,q) que mejor puede describir la serie. 3. Estimar los parámetros de los polinomios de rezagos (L) y (L) 4. Hacer un diagnóstico para con…rmar si es que el modelo utilizado en la estimación es consistente con la data. Este proceso de diagnóstico se basa fundamentalmente en las propiedades del residuo obtenido luego de estimar el modelo ARMA(p,q). Si el modelo asumido es el correcto, entoces este residuo debería no presentar correlación serial de ningún order, ser homocedástico, e idealmente debería tener media cero y distribuirse normalmente. De no cumplirse algunas de estas propiedades del residuo, se regresa al paso (2) y se repite el proceso. Adicionalmente, se puede utilizar criterios de información para escoger el modelo más adecuado. Este será el modelo que minimice el mayor número de criterios de información disponibles. Solamente con la …nalidad de precisar, reproducimos a continuación las fórmulas de los criterios de información más comúnmente usados. El criterio de información de Akaike ( AIC) AIC =
2 T L + 2 T k
El criterio de Schwarz (SC) SC =
2 T L + k log(T ) T
Y el criterio de Hanna-Quinn ( HQ) HQ =
2 T L + 2k log(log(T )) T
4.1 LA METODOLOGíA DE BOX Y JENKINS
27
En donde, T representa el tamaño de la serie de tiempo, L el logartimo de la función de verosimilitud, mayor este valor, mayor la probabilidad de que los datos prevengan del modelo planteado, y …nalmente, k que representa el número de parámetros estimados. Note que todos estos criterios son decrecientes en el valor de la función de verosimilitud, normalizados por el tamaño muestral, y crecientes en el número de parámetros estimados. En ese sentido, los criterios de información se basan en la misma idea que el R2 ajustado, es decir, medir el mejoramiento en el grado de ajuste pero eliminando el efecto que genera la inclusión de parámetros adicionales.
Capítulo 5
Estimación En las secciones anteriores aprendimos a derivar las propiedades, principalmente sus momentos de primer y segundo orden, de un conjunto de modelos de series de tiempo, asumiendo que los parámetros de los modelos analizados eran conocidos. En este capítulo, por el contrario, vamos a aprender como se pueden estimar estos parámetros de estos modelos a partir de los momentos muestrales, aquellos obtenidos con la data, de las series de tiempo que queremos analizar. El principio de estimación que vamos a utilizar es el de máxima verosimilitud. Un estimador de MV es aquel que maximiza la probabilidad conjunta de que la muestra observada haya sido generada por el modelo propuesto. Supongamos que tenemos una muestra de datos de tamaño T ; fy1 ; y2 ;::yT g, y proponemos que esta muestra ha sido generada por un modelo AR(1). yt = 1 yt1 + "t
donde, "t N (0; 2" ) El estimador de MV de 1 y 2" es aquel que maximiza la probabilidad conjunta L , a la que denominaremos función de verosimilitud, dada por: L = f Y T Y T
1
::Y 1 (y1 ;y2 ;::yT ;1 ; 2" )
Para encontar el estimador de MV, se requieren dos pasos. El primero consiste en encontrar la función de verosimilitud, y el segundo en hallar el valor de los parámetros que hace máxima esta función. Note que la función de verosimilitud es función tanto de los parámetros como de los datos muestrales, y por ello èsta deber estar bien de…nida respecto a estos dos argumentos. En particular, vamos a asumir que el espacio de parámetros sobre el cual se va optimizar es un espacio compacto, y que la función de 29
30
CAPíTULO 5 ESTIMACIÓN
verosimilitud es continua respecto a los parámetros. Estas dos condiciones garantizan la existencia de un máximo único. La función de verosimilitud puede escribirse de manera más sencilla en algunos casos particulares. Por ejemplo, bajo el supuesto que yt son variables independientes, que no es el caso del modelo AR(1) anterior, la función de verosimilitud puede escribirse como el producto de las funciones de densidad de cada variable aleatoria y t T
L(Y; ) =
Y
Lt (yt ; )
k=0
Aún en el caso que las variables no sean independientes, es posible, escribir la función de verosimilitud de una forma más compacta, utilizando las contribuciones marginales que vienen de cada observación . De…niendo como Lt (yt =yt1 ; yt2 ::;y1 ; ) la distribución de densidad condicional de yt , entonces la función de verosimilitud puede escribirse como: T
L(Y; ) =
Y
t=0
Lt (yt =yt1 ; yt2 ::;y1 ; )
Esta última forma de escribir la función de verosimilitud es la que se utiliza para estimar modelos de series de tiempo, puesto que en este caso el supuesto de independencia no se cumple. Es usual que se utilice para el proceso de estimación no directamente la función de verosimilitud sino una transformación monotónica de la misma, el logaritmo de la función de verosimilitud. T
log L(Y; ) =
X t
log Lt (yt =yt1 ; yt2 ::;y1 ; )
En la siguiete sección vamos a mostrar como calcular la función de verosimilitud para modelos ARMA(p,q), mientras en la subsiguiente, se hace una revisión rápida de los principales algoritmos de optimización disponibles para encontrar los estimadores de MV.
5.1. La función de verosimilitud El primer caso a revisar es el de un AR(1) Y t = c + Y t1 + "t
31
5.1 LA FUNCIÓN DE VEROSIMILITUD
donde, "t N (0; 2" ). De la sección anterior sabemos que los momentos no condicionales de Y t son: E (Y t ) = = 1c y E (Y t )2 = 1 " Sabemos de la de…nición de la función de verosimilitud que esta se puede obtener a partir de la función condicional de los T-1 variables y de la distribución marginal de Y 1. Puesto que "t se distribuye N (0; 2" ) tiene que ser el caso que Y 1 N (; 1" ), por lo que el logaritmo de su función de densidad esta dado por: 2
2
2
2
log L1 (y1 ; ) =
1 log(2) 2
1 2" log 2 1 2
1 (y1 )2 2" 2 2
1
Luego tenemos que calcular la función de densidad de Y 2 condicional en el valor de Y 1, dado el supuesto de normalidad, tenemos: L1 (y2 =y1 ; )
N (c + Y 1; 2" )
por lo que su función de verosimilitud estará dada por:
log L1 (y2 =y1 ; ) =
1 log(2) 2
1 log 2" 2
c Y 1)2
1 (y2 2
2"
Repitiendo este procedimiento hasta yT , y sumando todas las funciones de densidad relevantes, se obtiene la función de verosimilitud que tiene la siguiente forma: log L(Y; ) =
2
2
12 log(2) 12 log 1 12 (y ) "
1
2
2 "
12
T 1 (T 1) 2
log(2) (T 1) log 2 1 " 2 2
X t
(yt c Y t 2"
1
)2
Note que cuando T es grande, la contribución de la densidad marginal de la primera observación tiene es proporcionalmente menor. Esta es la razón por la que cuando el tamaño muestral es su…cientemente grade se puede prescindir de la primera observación, esto es, se asume esta como dada. En este caso, maximizar la función de verosimilitud es equivalente a maximizar, T 1
(y c Y )2 t t1
X 1 2
t
2"
32
CAPíTULO 5 ESTIMACIÓN
que no es otra cosa el negativo de la suma de cuadrados residuales. Entonces, el estimador de máxima verosimilitud obtenido utilizando la función de verosimilitud conditional es equivalente al estimador de mínimos cuadrados ordinarios. Puesto que computar la función de verosimilitud exacta para un AR(p) es un proceso relativamente complejo, para este tipo de procesos vamos a considerar la función de verosimilitud condicional 1 T 1
log L(Y; ) =
(T 1) 2
log(2) (T 1) log 2 1 " 2 2
X
(yt c1 Y t
1
t
Y t 2 2
2
::p Y t
"
Igualmente, para el caso de un modelo MA(1) vamos a considerar la función de verosimilitud condicional. Dado el siguiente proceso MA(1) Y t = + "t + "t1
que
donde, "t N (0; 2" ). Si el valor de " t1 es conocido, entonces, tenemos
N ( + "t1; 2" )
Y t ="t1
supongamos que conocemos que "0 = 0 entonces podemos calcular, "1 = Y 1
utilizando "1 podemos calcular "2 = Y 2 "1 y así sucesivamente podemos construir toda una secuencia de "t "t = Y t
"t1
Entonces, el log de la función de verosimilitud condicional para este proceso puede escribirse de forma muy sencilla como: T 1 2 1) (T 1) 1 "t 2 log L(Y; ) = log(2) log " 2 2 2 2
(T
X t
"
Igualmente, este procedimiento puede extenderse para el caso de un proceso MA(1), para calcular la función de verosimilitud condicional yt = + "t + 1 "t1 + 2 "t2 ::q "tq 1
En el capítulo de Filtros de Kalman se discutirá una forma alternativa que facilita el cálculo de la función de verosimilitud exacta
2
p
)
5.2 ALGORITMOS NUMÉRICOS DE OPTIMIZACIÓN 3
33
asumimos que: "o = "1 = :::"q+1 = 0
do:
Entonces recursivamente, se puede obtener una secuencia de " t utilizan"t
yt 1"t1 2"t2:: q "tq
Finalmente, tenemos que:
T 1 2 1) (T 1) 1 "t 2 log L(Y; ) = log(2) log " 2 2 2 2
X
(T
t
"
5.2. Algoritmos numéricos de optimización2 En esta sección vamos a evaluar de manera suscita una serie de algoritmos numéricos de optimización que luego nos van a permitir hacer estimación mediante el método de máxima verosimilitud. Vamos a revisar dos tipos de algoritmos, aquellos que usan derivadas y aquellos que no. El problema típico de optimización puede representarse como : m ax f (x) xX
donde f es la función objetivo y X es el conjunto factible, y x representa el máximo si este existe. El teorema de Wiertrass nos ofrece las condiciones bajo las cuales existe una solución para este problema. Si f es una función contínua y X es no vacío, cerrado y acotado, entonces f tiene un máximo en X . Una condición necesaria y su…ciente para encontrar un máximo global es que la función f sea estrictamente cóncava. Todos los métodos númericos de optimización se basan en evaluaciones contínuas de la función objetivo en diferentes puntos del conjunto factible, la diferencia entre ellos radica fundamentalmente en como se acercan a la solución dentro del conjunto factible, aquellos más e¢cientes lo hacen más rápido utilizando información sobre el gradiente, o la Hesiana de la función. Sin embargo en muchos casos, calcular el gradiente puede ser una tarea complicada por lo que algoritmos simples que no utilicen mayor información pueden resultar más convenientes. 2
Este capítulo se basa en el Libro .A pplied Computational Economics and Finance"de Mario J. Miranda y Paul L. Fackler
34
CAPíTULO 5 ESTIMACIÓN
Todo estos métodos iteran la búsqueda de la solución hasta que se alcance algún criterio de convergencia. Puesto que no existe una solución analítica no es posible determinar la solución exacta, se tiene que aceptar un grado de aproximación. Este esta dado por el criterio de convergencia, que usualmente consiste en un cambio arbitrariamente pequeño en la función objetivo o en la solución de la misma. Que tan pequeño sea este cambio determina la precisión de la solución.
5.2.1. Métodos libres de derivadas Estos médos utilizan evaluaciones sucesivas de la función en una partición …na X , para luego escoger el punto que hace f (x ) máximo. El método más sencillo es el grid search, que no realiza ningun tipo de ajuste al procedimiento de búsqueda. Simplemente se evalúa la función objetivo en una partición muy …na de conjunto factible, luego se escoge el valor de x que hace la función f máxima. Ejercicio: Escriba un código en Matlab que permita determinar el máximo de la función f = 1 + x 0;5x2 utilizando el método de grid search.
5.2.2. Métodos basados en derivadas La principal ventaja de estos métodos es su mayor rápidez para encontrar la solución. Puesto que usan información de la primera y segunda derivada permite a…nar la dirección de búsqueda de la solución. En general todos estos métodos utilizan la siguiente regla para acercarse a la solución:
b b
s+1 = s + s As gs
b
donde As es un matriz que depende de los parámetros estimados, s , gs es el gradiente de la función evaluada también en s , y s es un parámetro que calibra la velocidad de ajuste en la dirección del gradiente. Para poder implementar métodos de este tipo, usualmente se requiere un punto inicial 0 . La elección del punto de partida para la búsqueda es muy importante para garantizar que se pueda alcanzar la solución de manera rápida. Más aún el no escoger apropiadamente los valores iniciales puede hacer que los métodos de búsqueda no funcionen. Cuando la función no es globalmente cóncava, es siempre buena idea, utilizar varios puntos iniciales, para chequear si el máximo obtenido es un máximo local. Además se requiere, de…nir un criterio de convergencia, y s y A s . Dependiendo como se escojan estos parámetros tenemos diferentes métodos de optimización:
b
b
5.2 ALGORITMOS NUMÉRICOS DE OPTIMIZACIÓN 4
35
Diferentes Métodos de Optimización Método
s
As
Newton Raphson
1
BHHH
H s = @ @f
1
Ascenso rápido
0
gs gs gsH s gs 0
2
2
H HH =
@f @
@f @
0
I q
Puesto que para implementar cualquiera de estos algoritmos de optimización se requiere calcular derivadas, el procedimiento puede complicarse de manera importante a menos que las derivadas se puedan calcular numéricamente.
Capítulo 6
Aplicaciones de Laboratorio 6.1. Identi…cación de los procesos ARIMA 6.1.1. AR(1) Viene dado por donde jj < 1,
yt = + yt1 + "t E (yt ) = y =
(6.1)
1
2 E y t = 2y = 0 = k =
2" 1 2
k = k ; k > 0 0
8
Los procesos AR(1) se reconocen por una FAS in…nita y una FAP que se desvanece tras el primer rezago. Es necesario especi…car un término constante a siempre que los datos tengan media distinta de cero
. Función de autocorrelacion simple y parcial (resp.) para un AR(1) cuando 0 < < 1 .
37
38
CAPíTULO 6 APLICACIONES DE LABORATORIO
6.1.2. MA(1) Viene dado por yt = y + "t
, jj < 1
"t1
E (yt ) = y ;
2 = 2 = = 1 + 2 2 y 0 "
E y t
k = k = 0
8< :
; k = 1+2
1
0; k > 1
Los procesos MA(1) se reconocen por una FAP in…nita y una FAS que se desvanece tras el primer rezago. Es necesario especi…car un término constante en los datos tengan media distinta de cero.
Función de autocorrelacion simple y parcial (resp.) para un MA(1) cuando 0 < < 1 .
6.1 IDENTIFICACIÓN DE LOS PROCESOS ARIMA
39
Función de autocorrelacion simple y parcial (resp.) para un MA(1) cuando 1 < < 0 .
6.1.3. AR(2) Viene dado por yt = + '1 yt1 + '2 yt2 + "t ; '1 + '2 < 1 , '2
'1 < 1, j'2j < 1
E (yt ) = y =
'1
'2 ; k = k = ' 1 k1 + '2 k2 ; k 1; 1
0
donde: 1 =
'1 '21 ; 2 = ' 2 + 1 '2 1 '2
Los procesos AR(2) se reconocen por una FAS in…nita y una FAP que se desvanece tras el segundo rezago. Es necesario especi…car un término constante en tanto los datos tengan media distinta de cero.
Función de autocorrelacion simple y parcial (resp.) para un AR(2) cuando ' 1 > 0 , ' 2 > 0 .
40
CAPíTULO 6 APLICACIONES DE LABORATORIO
Figura 6.1: Función de autocorrelacion simple y parcial (resp.) para un AR(2) cuando '1 < 0, ' 2 > 0 .
Función de autocorrelacion simple y parcial (resp.) para un AR(2) cuando ' 1 > 0 , ' 2 < 0 .
Función de autocorrelacion simple y parcial (resp.) para un AR(2) cuando ' 1 < 0 , ' 2 < 0 .
6.1.4. MA(2) Viene dado por yt = y + "t
1"t1 2"t2
6.1 IDENTIFICACIÓN DE LOS PROCESOS ARIMA
41
Figura 6.2: Función de autocorrelacion simple y para un MA(2) cuando 1 < 0 , 2 > 0 y 1 > 0, 2 < 0 (resp.).
1 + 2 < 1; 2
2
E y t
1 < 1; j2j < 1
E (yt ) = y ;
= 2y = 0 = 1 + 21 + 22 2" ;
k =
k = 0
8>< >:
(1 ) ; k = 1 1+ + ; k = 2 1
2 2 2
2 1
2
1+21 +22
0; k > 2
Los procesos MA(2) se reconocen por una FAP in…nita y una FAS que se desvanece tras el primer rezago. Es necesario especi…car un término constante en tanto los datos tengan media distinta de cero.
Función de autocorrelacion simple y para un MA(2) cuando 1 > 0, 2 > 0 y 1 < 0, 2 < 0 (resp.).
6.1.5. Comandos usuales en Eviews Los comandos que suelen utilizarse para este tipo de evaluaciones son los siguientes:
42
CAPíTULO 6 APLICACIONES DE LABORATORIO
Regresión Mínimo Cuadrática (MCO): LS Se indica la variable endógena seguidamente de las explicativas: Ej: LS y c x1 : y t = + x 1t + "t . o Para incluir un rezago como explicativa indicar AR(p): Ej: LS y c AR(1) : yt = + '1 yt1 + "t . o Para incluir un rezago del error como explicativa indicar MA(q): Ej: LS y c MA(1) : y t = y + "t 1 "t1 . o Para incluir una tendencia determinística como explicativa indicar
@trend: Ej: LS y c @trend : yt = + t + "t . o Para incluir una raíz estacional como explicativa indicar @seas(i): Ej: LS y c @seas(1) : y t = + '1 yt4 + "t Visualización de Grá…cos y Tablas: show Se indica el nombre de la serie seguidamente del tipo de grá…co a mostrar: Show x1.line : Visualización del grá…co lineal de la serie x1. Show x1.ident(24) : Visualización del correlograma de la serie x1 con 24 rezagos. Show d(x1).ident(24) : Visualización del correlograma con 24 rezagos de la primera diferencia de la serie x1.
6.1.6. Apéndice: Simulación de procesos ARMA en MATLAB Esta sección está dedicada a reforzar los conceptos adquiridos hasta este momento en clase. Si bien es cierto el Eviews nos permite observar fácilmente las características principales de los procesos ARMA, no es del todo claro a la hora de tratar de entender cómo se construyen las cosas. El curso de Series de Tiempo está pensado para desarrollar conceptos teóricos, así como para llevarlos a la practica, apuntando a que no queden solamente en el papel. En suma, lo que se pretende es implementar en el computador las principales derivaciones analíticas vistas en clase, desde la generación de un proceso ARMA cualquiera hasta la construcción de las FAS, FAP y la Función de Impulso Respuesta. Para ello haremos uso del paquete MATLAB, que es escogido entre muchos otros por su versatilidad con el tabajo de matrices y operadores lógicos, así como por su alta potencialidad, lo cual permite trabajar con algoritmos de optimización a gran escala, que serán de mucha utilidad en el desarrollo del curso. Cabe señalar que en la mayoría de capítulos se incluirá un apartado similar, apuntando a reforzar lo anteriormente señalado.
6.1 IDENTIFICACIÓN DE LOS PROCESOS ARIMA
43
Las características principales de estos procesos se encuentran en la sección 1.2, por lo que en esta oportunidad no entraremos en dichos detalles.
AR(1) Viene dado por y t = + yt1 + "t : El código MATLAB para generar este proceso se presenta a continuación. Este archivo .m es de tipo función, lo cual implica que para que funcione tenemos que asignarle valores u objetos a los argumentos del mismo. ar1.m:
f unction [y FAS]=ar1(c,rho,T,sigma)
% c: Término constante. % rho: Parámetro de persistencia. % T: Número de observaciones del experimento. % sigma: Desviación estándar del término de error. y0=randn(1,1); % Especi…camos la condición inicial, proveniente de una distribución normal. y(1,1)=c+rho*y0(1,1)+sigma*randn(1,1);% Generamos la primera observación. for i=2:T y(i,1)=c+rho*y(i-1,1)+sigma*randn(1,1);% A partir de la primera generamos las T-1 restantes. end Mu=c/(1-rho);% valor esperado del proceso generado. g0=(sigma^2)/(1-rho^2); %Varianza del proceso generado. L=25; for j=1:L g(j,1)=(rho^j/(1-rho^2))*(sigma^2); % Función de Autocovarianza FAS(j,1)=g(j,1)/g0; % Función de Autocorrelación end
Una vez creada la función que genera el proceso, realizamos un experimento para distintos valores de persistencia. La función anterior calcula también la función de autocorrelación simple (FAS) y para el cálculo de la Función de impulso respuesta utilizamos la función …lter, que será detallada con el desarrollo del curso. Este nuevo archivo no es función, sino un procedimiento. Es decir, es un conjunto de comando que se realizan en una secuencia determinada, pero que no dependen de ningún argumento para poder funcionar. Ello no quita la posibilidad de que este tipo de archivos
44
CAPíTULO 6 APLICACIONES DE LABORATORIO
pueda incluir el llamado de cualquier otra función generada anteriormente. Tomando ventaja de ello tenemos: exp_ar1.m: % Experimentos con AR(1) T=100; c=0; sigma=1; m=zeros(25,1); m(1,1)=1; b=1; rho=0; [y(:,1) FAS(:,1)]=ar1(c,rho,T,sigma); a=[1 -rho]; f(:,1)=…lter(b,a,m); rho=0.5; [y(:,2) FAS(:,2)]=ar1(c,rho,T,sigma); a=[1 -rho]; f(:,2)=…lter(b,a,m); rho=0.9; [y(:,3) FAS(:,3)]=ar1(c,rho,T,sigma); a=[1 -rho]; f(:,3)=…lter(b,a,m); …gure(1) subplot(3,1,1);plot(y(:,1),’b’);xlabel(’White Noise (rho=0)’); subplot(3,1,2);plot(y(:,2),’b’);xlabel(’AR(1) (rho=0.5)’); subplot(3,1,3);plot(y(:,3),’b’);xlabel(’AR(1) (rho=0.9)’); …gure(2) subplot(3,1,1);BAR(FAS(:,1));xlabel(’FAS-White Noise (rho=0)’); subplot(3,1,2);BAR(FAS(:,2));xlabel(’FAS-AR(1) (rho=0.5)’); subplot(3,1,3);BAR(FAS(:,3));xlabel(’FAS-AR(1) (rho=0.9)’); …gure(3) subplot(3,1,1);BAR(f(:,1));xlabel(’Impulse-Response - White Noise (rho=0)’); subplot(3,1,2);BAR(f(:,2));xlabel(’Impulse-Response - AR(1) (rho=0.5)’); subplot(3,1,3);BAR(f(:,3));xlabel(’Impulse-Response - AR(1) (rho=0.9)’);
6.1 IDENTIFICACIÓN DE LOS PROCESOS ARIMA
6.1.7. AR(2) Viene dado por yt = + '1 yt1 + '2 yt2 + "t ; '1 + '2 < 1, ' 2
'1 < 1, j'2j < 1
45
46
CAPíTULO 6 APLICACIONES DE LABORATORIO
El procedimiento es similar al anterior, con la salvedad de que hay que especi…car dos condiciones iniciales en vez de una ar2.m function [y FAS]=ar2(c,rho1,rho2,T,sigma) y0=randn(2,1); y(1,1)=c+rho1*y0(2,1)+rho2*y0(1,1)+sigma*randn(1,1); y(2,1)=c+rho1*y(1,1)+rho2*y0(2,1)+sigma*randn(1,1); for i=3:T y(i,1)=c+rho1*y(i-1,1)+rho2*y(i-2,1)+sigma*randn(1,1); end Mu=c/(1-rho1-rho2); L=25; g0=(sigma^2)/(1-rho1^2-rho2^2); g(1,1)=(rho1*g0)/(1-rho2); g(2,1)=rho1*g(1,1)+rho2*g0; FAS(1,1)=g(1,1)/g0; FAS(2,1)=g(2,1)/g0; for j=3:L g(j,1)=rho1*g(j-1,1)+rho2*g(j-2,1); FAS(j,1)=g(j,1)/g0; end Un experimento con distintos valores es similar al anterior, haciendo pequeños cambios en el código.
6.1.8. MA(1) Viene dado por yt = y + "t "t1 , < 1 En este caso la función de autocorrelación es …nita, tal como vimos en la parte de identi…cación. ma1.m function [y FAS]=ma1(c,theta,T,sigma) e=sigma*randn(T,1); y(1,1)=randn(1,1); for i=2:T
y(i,1)=c+e(i,1)-theta*e(i-1,1);
jj
6.1 IDENTIFICACIÓN DE LOS PROCESOS ARIMA
47
end Mu=c; g0=(1+theta^2)*sigma^2; g(1,1)=-(theta/(1+theta^2))*sigma^2; FAS(1,1)=g(1,1)/g0; L=25; for j=2:L g(j,1)=0; FAS(j,1)=g(j,1)/g0; end
De igual forma, un experimento con distintos valores es similar al anterior, haciendo pequeños cambios en el código.
Proceso ARMA(p,q) generalizado Todos los proceso anteriores, así como el cálculo de sus respectivas FAS y Funciones de Impulso Respuesta pueden ser resumidos en un solo código, tal como sigue: armapq.m function [y FAS]=armapq(A,B,T) e=randn(T,1); y=…lter(B,A,e); L=25; FAS=autocorr(y,L); end
Del mismo modo, realizamos el experimento para diversos modelos: exp_armapq.m % Experimentos con ARMA(p,q) T=100; m=zeros(25,1); m(1,1)=1; %ARMA(1,1) A=[1 -0.7]; B=[1 0.5]; [y(:,1) FAS(:,1)]=armapq(A,B,T); f(:,1)=…lter(B,A,m); %ARMA(2,1)
48
CAPíTULO 6 APLICACIONES DE LABORATORIO A=[1 -0.7 0.2]; B=[1 0.5]; [y(:,2) FAS(:,2)]=armapq(A,B,T); f(:,2)=…lter(B,A,m); %ARMA(1,2) A=[1 -0.7]; B=[1 0.5 -0.2]; [y(:,3) FAS(:,3)]=armapq(A,B,T); f(:,3)=…lter(B,A,m); …gure(1) subplot(3,1,1);plot(y(:,1),’b’);xlabel(’ARMA(1,1) (rho1=0.7;theta1=-0.5)’);
subplot(3,1,2);plot(y(:,2),’b’);xlabel(’ARMA(2,1) (rho1=0.7;rho2=-0.2;theta1=0.5)’); subplot(3,1,3);plot(y(:,3),’b’);xlabel(’ARMA(1,2) (rho1=0.7;theta1=-0.5;theta2=0.2)’); …gure(2) subplot(3,1,1);BAR(FAS(:,1));xlabel(’FAS-ARMA(1,1) (rho1=0.7;theta1=-0.5)’); subplot(3,1,2);BAR(FAS(:,2));xlabel(’FAS-ARMA(2,1) (rho1=0.7;rho2=-0.2;theta1=0.5)’); subplot(3,1,3);BAR(FAS(:,3));xlabel(’FAS-ARMA(1,2) (rho1=0.7;theta1=-0.5;theta2=0.2)’); …gure(3) subplot(3,1,1);BAR(f(:,1));xlabel(’Impulse-Response - ARMA(1,2) (rho1=0.7;theta1=0.5;theta2=0.2)’); subplot(3,1,2);BAR(f(:,2));xlabel(’Impulse-Response - ARMA(1,2) (rho1=0.7;theta1=0.5;theta2=0.2)’); subplot(3,1,3);BAR(f(:,3));xlabel(’Impulse-Response - ARMA(1,2) (rho1=0.7;theta1=0.5;theta2=0.2)’);
6.1.9. Asignación Se presenta a continuación el archivo Clase1.wf1 que contiene 12 series de tiempo de 500 observaciones cada una.
6.2 PREDICCIÓN Y EXPERIMENTOS COMPUTACIONALES 49
Se le pide identi…car el tipo de proceso ARIMA(p,d,q) que sigue cada serie, utilizando los criterios de discriminación de modelos estudiados en clase. Se le pide: En grupos de 2 personas seleccione 6 de las 12 series que se encuentran en el archivo Clase1.wf1 e identi…que el tipo de proceso ARIMA que sigue cada una de ellas. Se le aconseja grabar una ecuación por cada serie analizada. En el objeto de texto (_datos) escriba el nombre de cada integrante del grupo, las 6 series seleccionadas y el proceso que sigue cada una de ellas. Modi…que el nombre del work…le asignando Clase1_I1_I2.wf1, siendo In el primer apellido del integrante n del grupo.
6.2. Predicción y Experimentos Computacionales 6.2.1. Introducción En este capítulo implementaremos los conceptos desarrollados en el anterior básicamente en MATLAB, sobretodo tratando de llegar al tema de evaluacion de resultados y calidad de los mismos. El desarrollo computacional de estos aspectos es un complemento necesario para el aprendizaje, en la medida en que permite esclarecer los criterios numéricos desarrollados en forma analítica. Cabe señalar que la versión de MATLAB utilizada es la 7.0 R14, que contiene el GARCH toolbox, necesario para varias aplicaciones que serán desarrolladas líneas abajo. Al …nal del capítulo se ofrecerá un apartado sobre predicción de series temporales utilizando Eviews 5.1.
50
CAPíTULO 6 APLICACIONES DE LABORATORIO
6.2.2. Lineamientos básicos Recordando la descomposición de Wold: yt = c +
1
X k=0
k "tk
(6.2)
yt = c + "t + 1 "t1 + 2 "t2 ::::
Supongamos que y es un proceso ARMA ( p; q ) generado con el código armapq.m . Luego la función garchma es utilizada para convertir dicho proceso en un M A (1) de la siguiente forma:
phi = garchma(AR,MA,k);
Donde phi es el vector de k componentes. mostrado en la ecuación anterior. Si k no es especi…cado, por defecto MATLAB tomará k = 10. Si bien es cierto la relación teorica demanda un término in…nito (1), numéricamente es necesario …jar un criterio de convergencia tal como k , de modo que a partir de ese valor los rezagos no añadan valor a la serie. Cabe señalar que (1) = 1, tomando como referencia al término que premultiplica a la perturbación en el período t .AR es el vector de dimensión p que contiene a los parámetros de tipo autorregresivo del modelo. MA es análogo a AR, de dimensión q y que contiene a los parámetros de medias móviles. Es importante mencionar que los vectores AR y MA deben cumplir con las condiciones de estacionariedad e invertibilidad respectvamente. Contando con los parámetros de la descomposición de wold, es posible replicar las fórmulas del capítulo anterior: Varianza del error de predicción:
var(yt ) =
2
1
X
2k
k=0
I=ones(k,1); Var_y=(sigma^2)*(I*phi);
Error cuadrático medio:
MSE (t+h=t ) = 2 1 + 21 + :: + 2h1
6.2 PREDICCIÓN Y EXPERIMENTOS COMPUTACIONALES 51 phi2=phi.^2; phi2=phi2(1:k-1); I=ones(k-1,1); MSE=(sigma^2)*(phi2*I);
Bandas de con…anza: :yt+h=I t
q
1;96
I=ones(T,1); yhtb1=yht + 1.96*sqrt(MSE); yhtb2=yht - 1.96*sqrt(MSE);
2 1 + 21 + :: + 2h1
6.2.3. Predicción y bandas de con…anza en un ARMA ( p; q ) El ejemplo que sigue es la simulación de un proceso ARMA ( p; q ) y el análisis de predicción del mismo: % Predicción con Modelos ARMA(p,q) estacionarios % % Fernando Pérez Forero % 04/12/05 clear % Generamos un proceso ARMA(p,q) c=0;% Intercepto A=[1 -0.6 0.2]; % Componente AR B=[1 0.5];% Componente MA T=1000;% Número de Obs. sigma=1;% Desv. Estándar del término de perturbación [y e]=armapq(A,B,T,sigma); y=c+y; % Generamos el vector de parámetros de la Descomposición de Wold para este % proceso. n=length(A); AR=A(2:n); AR=-1*AR; MA=B; k=15; phi=garchma(AR,MA,k); % Vector de parámetros de la R. Wold. I=ones(k,1); Var_y=(sigma^2)*(I*phi);% Varianza del proceso y.
52
CAPíTULO 6 APLICACIONES DE LABORATORIO % Predicción h períodos hacia adelante. h=2; beta=phi(h:k); s=length(beta); for i=T:-1:h+s ep=e(i-s+1:i); yht(i,1)=c+beta*ep; % y(t+hjt) end
% Error Cuadrático Medio (MSE) phi2=phi.^2; phi2=phi2(1:k-1); I=ones(k-1,1); MSE=(sigma^2)*(phi2*I); % Bandas de Con…anza I=ones(T,1); yhtb1=yht + 1.96*sqrt(MSE); yhtb2=yht - 1.96*sqrt(MSE); % Grá…co …gure(1) plot(yht,’g’);title(’ARMA(2,1) - Actual vs. Forecast’); hold on plot(y,’b’) hold o¤ …gure(2) plot(yhtb1,’r’);title(’Predicción ARMA(2,1) - Bandas al 95 %’); hold on plot(yht,’b’) hold on plot(yhtb2,’r’) hold o¤
6.2 PREDICCIÓN Y EXPERIMENTOS COMPUTACIONALES 53
6.2.4. Encontrando el mejor predictor lineal en un AR( p) Sea un proceso AR( p) como sigue: (L) (yt
donde: (L) = 1
c) = t
1L 2L2:::: pL p
El algoritmo para la predicción de un modelo AR(p) es esencialmente el mismo que un AR(1), puesto que un modelo AR(p) se puede escribir como un modelo AR(1) utilizando la representación de estado espacio.
0B BB @
yt
1C CC A
c yt1 c : yt p+1
c
=
0B B@
1 2 :: p 1 0 :: 0 0 :: 1 0
1C 0B CA BB@
yt1
1C CC A
c yt2 c : yt p
c
+
0B BB @
t 0 : 0
1C CC A
54
CAPíTULO 6 APLICACIONES DE LABORATORIO o en términos de notación matricial, se tiene que: t = F t1 + wt var(wt ) = w
% Representación de Estado-Espacio F=[AR]; F=[F;eye(p-1),zeros(p-1,1)]; for i=3:T E(:,i)=[y(i,1);y(i-1,1);y(i-2,1)]; w(:,i)=[e(i,1);zeros(p-1,1)]; end
Comenzando en el periodo t, y iterando hacia adelante h periodos, tenemos:
t+h = F h t + wt+h + F wt+h1 :::::F h1 wt+1
% h=1 for i=4:T Eh1(:,i)=(F^1)*E(:,i-1)+w(:,i)+(F^1)*w(:,i-1); end % h=2 for i=5:T Eh2(:,i)=(F^2)*E(:,i-2)+w(:,i)+(F^1)*w(:,i-1)+(F^2)*w(:,i-2); end
Entonces el mejor predictor lineal basado en información hasta I t t+1=t = F t t+2=t = F 2 t t+h=t = F h t
% Hallando el mejor predictor lineal % h=1
6.2 PREDICCIÓN Y EXPERIMENTOS COMPUTACIONALES 55 for i=4:T Eh1t(:,i)=F*E(:,i-1); end % h=2 for i=5:T Eh2t(:,i)=(F^2)*E(:,i-2); end
El error de predicción asociado están dados por: wt+1=t = t+1 wt+2=t = t+2 wt+h=t = t+h
t+1=t = wt+1
t+2=t = wt+2 + F wt+1
t+h=t = wt+h + F wt+h1 + :: + F hwt+1
% Error de predicción err1=Eh1- Eh1t; err2=Eh2- Eh2t;
y el correspondiente matrices de varianza-covarianza del error de predicción. var(wt+1=t ) = w var(wt+2=t ) = w + F w F 0 h 1
var(wt+h=t ) =
X
F j w F j 0
j=0
% Varianza del error de predicción varw=w*w’; H=20; varwh=varw; for i=1:H for h=1:i varwh(:,:,i)=varw; varwh(:,:,i)=varwh(:,:,i)+(F^(h-1))*varw*(F^(h-1))’; end end
56
CAPíTULO 6 APLICACIONES DE LABORATORIO Notese que: var(wt+h=t ) = w + Fvar(wt+h=t )F 0
Para una mejor compresion del ejemplo presentamos el ejercicio completo a continuación: % Predicción de un AR(P) % % Fernando Pérez Forero % 11/12/05 % clear % Generamos un Proceso AR(P) P=3 c=0;% Intercepto A=[1 -0.6 0.2 -0.5];% Componente AR B=[1]; % Componente MA T=1000; % Número de Obs. sigma=1; % Desv. Estándar del término de perturbación [y e]=armapq(A,B,T,sigma); y=c+y; % Capturamos los parametros asociados a los rezagos n=length(A); AR=A(2:n); AR=-1*AR; p=n-1; % Representación de Estado-Espacio F=[AR]; F=[F;eye(p-1),zeros(p-1,1)]; for i=3:T E(:,i)=[y(i,1);y(i-1,1);y(i-2,1)]; w(:,i)=[e(i,1);zeros(p-1,1)]; end % Predicción h períodos hacia adelante % h=1 for i=4:T Eh1(:,i)=(F^1)*E(:,i-1)+w(:,i)+(F^1)*w(:,i-1);
6.2 PREDICCIÓN Y EXPERIMENTOS COMPUTACIONALES 57 end % h=2 for i=5:T Eh2(:,i)=(F^2)*E(:,i-2)+w(:,i)+(F^1)*w(:,i-1)+(F^2)*w(:,i-2); end % Hallando el mejor predictor lineal % h=1 for i=4:T Eh1t(:,i)=F*E(:,i-1); end % h=2 for i=5:T Eh2t(:,i)=(F^2)*E(:,i-2); end % Error de predicción err1=Eh1- Eh1t; err2=Eh2- Eh2t; % Varianza del error de predicción varw=w*w’; H=20; varwh=varw; for i=1:H for h=1:i varwh(:,:,i)=varw; varwh(:,:,i)=varwh(:,:,i)+(F^(h-1))*varw*(F^(h-1))’; end end
6.2.5. Predicción en Eviews El programa Eviews cuenta con una sentencia que ejecuta el procedimiento de predecir los valores sobre la muestra. Esta vista nos muestra los indicadores básicos de predicción, tales como la U de Theil, los sesgos de
58
CAPíTULO 6 APLICACIONES DE LABORATORIO
varianza y covarianza. Estos resultados se pueden observar dentro de un objeto equation, dado que resulta de una regresión mínimo cuadrática. Antes de realizar dicha tarea, debemos repasar lo hecho en clase sobre la forma en que realmente predecimos los valores de un proceso ARIMA(p,d,q). Para ello necesitamos recordar conceptos como predictores lineales, conjuntos de información disponible, error cuadrático medio, etc. Comenzaremos con un ejercicio similar al de la clase anterior, donde generamos la serie que vamos a luego predecir. La serie que vamos a analizar sigue el proceso: yt = y + yt1 + "t
donde , t N (0; 2 ) En base a esto, la mejor predicción lineal es: yt+hjt = y + h (yt
)
Donde el error de predicción correspondiente es: "t+hjt = y t+h
yt+hjt
Y la varianza:
V ar "t+hjt
1 2h = (1 2 ) 2
Donde debe cumplirse que:
l{m yt+hjt = y = E (yt )
h
!1
l{m V ar "t+hjt
h
!1
2 = = V ar (yt ) (1 2 )
Trabajando en Eviews primero generamos los parámetros del programa: ’Predicción de un AR(1) h pasos hacia adelante ’Fernando Pérez Forero ’01/07/05 ’Número de observaciones !Obs=100 ’Parámetros !mu=0.5 !rho=0.8
6.2 PREDICCIÓN Y EXPERIMENTOS COMPUTACIONALES 59 !sigma2=1 ’H-step ahead !h=7
Luego, creamos un work…le y generamos la serie AR(1): work…le u !Obs series AR1=0 series z=nrnd smpl @…rst+1 @last series AR1=!mu+!rho*AR1(-1)+z smpl @all
Luego, creamos las series que almacenarán tanto los valores predichos como las respectivas varianzas: series AR1f=AR1 series var_e=0 smpl @last-!h+1 @last series AR1f=0 smpl @all
Con todo esto, comenzamos a construir los valores predichos para las h últimas observaciones: for !i=1 to !h AR1f(!Obs-!h+!i)=!mu+(!rho^!i)*(AR1(!Obs-!h)-!mu) var_e(!Obs-!h+!i)=!sigma2*(1-!rho^(2*!i))/(1-!rho^2) next
Computamos los errores de predicción, dado que contamos con los valores actuales: series e_forecast_…nal=AR1f-AR1 smpl @last-!h+1 @last
Construimos los grá…cos para comprobar las hipótesis de convergencia y comportamiento: freeze(e_forecast) e_forecast_…nal.line show e_forecast freeze(var_error) var_e.line show var_error freeze(ar1_f) ar1f.line show ar1_f smpl @all
60
CAPíTULO 6 APLICACIONES DE LABORATORIO
6.2.6. Análisis de indicadores La evaluación de las predicciones se hacen también dentro del objeto equation, oprimiendo el botón ‘Forecast’. En el resultado se observan los valores de la U de Theil, sesgo de media, varianza y covarianza, cuyas de…niciones se presentan a continuación: U de Theil:
v ut X v ut X v ut X T +n
^
yt
U =
yt
2
=n
t=T +1
T +n
^2
yt =n +
t=T +1
T +n
yt2 =n
t=T +1
Sesgo de Media: Cuando la media de la predicción dista mucho de la
de los datos corrientes. Debe ser cercano a cero, para ser considerada una predicción aceptable. Sesgo de Varianza: Cuando la variabilidad de la predicción no es la misma de la de los datos corrientes, teniendo problemas para explicar movimientos de la serie. Es preferible que este sesgo sea pequeño, de forma que las predicciones de corto plazo se encuentren bien alineadas con los datos reales, en la medida en que el modelo captura dicho movimiento.
6.3 ESTIMACIÓN POR MÁXIMA VEROSIMILITUD
61
Sesgo de Covarianza: Cuando la variabilidad de la predicción es bastante similar, pero no en el instante exacto en que indican los datos corrientes. Es preferible que el mayor sesgo de predicción se ubique en este segmento, dado que estamos capturando tanto tendencia como variabilidad, siendo la diferencia un sesgo de ocurrencia simultánea de eventos.
6.2.7. Ejercicios 1. En grupos de 2 personas, construya las bandas de con…anza del error de predicción al 95% de con…abilidad utilizando el programa que se ha construido en clase. Luego entregue el archivo al profesor. Pistas: El valor de la función de densidad de la distribución t es 1,96. La construcción de la banda es : yt+hjt =
1;96
2. Generalice el programa anterior a un proceso AR(2).
6.3. Estimación por Máxima Verosimilitud 6.3.1. Introducción Luego de revisar el tema de identi…cación de los procesos ARIMA, la tarea siguiente es la de estimar correctamente los parámetros de dicho modelo detectado. Para ello vamos a tener en cuenta la metodología desarrollada
62
CAPíTULO 6 APLICACIONES APLICACIONES DE LABORATORIO LABORATORIO
en clase, que es la de Máxima Verosimilitud. En estricto, la estimación por Máxima Verosimilitud ejecuta un procedimiento en el cual se debe asumir una distribución teórica a priori, en otras palabras, el estimador será de tipo paramétrico. Siendo este el escenario, buscaremos evaluar si la distribución asumida es la correcta mediante los tests adecuados, así como la e…ciencia de la estimación.
6.3.2. Algoritmo Algoritmoss de Optimizac Optimización ión y Máxim Máximaa Verosim Verosimilitud ilitud El procedimiento que vamos a replicar en este capítulo va va de la mano con el tema de la Optimizacion Matemática de funciones respecto a un número …nito de variables. Es decir, utilizaremos un algoritmo que se encargará de encontrar valores máximos para una función f de…nida en un dominio Rn. Dicho algortimo encontrará X encontrará un valor valor x 2 X tal tal que f (x (x ) f (x (x) para cualquier x 2 X . De…nimos el problema como: m ax ax f (x (x) x X
2
y llamamos a f la función objetivo a maximizar. Este tipo de optimización (de dimensión …nita) es frecuentemente frecuentemente utilizado en el análisis económico. económico. Por ejemplo, en los modelos neoclásicos estándar de la …rma y del consumidor se requiere maximizar tanto una función de bene…cios como una de utilidad respectivamente.Los modelos equilibrio competitivo son frecuentemente frecuentemente caracterizados caracterizados como problemas de optimización optimización donde existe un plani…cador social que se encargará de tomar las decisiones de todos los sectores en forma conjunta, de modo que maximice el bienestar de la sociedad. En general podemos encontrar diferentes algoritmos que se encargan de llevar a cabo este proceso, pero debemos tener en cuenta que ciertas características especiales pueden hacer variar un procedimiento con respecto a otro, tales como la aproximación de la matriz hessiana, el criterio de convergencia, así como el método para avanzar avanzar desde las condiciones iniciales hasta la respuesta …nal, de forma tal que sea un máximo global y no local. Tanto Eviews como Matlab cuentan con rutinas especializadas en estos tópicos, de manera que no tendremos mayores di…cultades a la hora de realizar los experimentos computacionales. computacionales. En econometría, la primera noción de optmizacion la tenemos al efectuar una estimación por Mínimos Cuadrados Ordinarios (MCO), dado que hallamos los parámetros para los cuáles la suma de errores al cuadrádo es mínima. Como se señaló en la introducción, este tipo de estimación es no paramétrica. El mismo procedimiento puede ser replicado si asumimos un
6.3 ESTIMACIÓN POR MÁXIMA VEROSIMILITUD
63
tipo de distribución particular, particular, que nos permita cuanti…car cuanti…car las densidades de probabilidad de cada realización del proceso. La obtención de los paramtros que maximizan el producto de densidades de probabilidad se denomina proceso de Máxima Verosimilitud. En estricto, la estimación por Máxima Verosimilitud es un método para obtener estimadores de los parámetros de un modelo, que tienen una distribución especí…ca. A diferencia de la estimación por MCO, que solamente toma en cuenta los dos primero momentos de la distribución, Máxima Verosimilitud incorpora toda la informacion de la distribución en el modelo, dado que trabaja con la distribución conjunta de las observaciones. El procedimiento es como sigue: Sea f una función de densidad de probabilidades, tal que f yt jt1 es la densidad condicional del proceso yt en el instante t, sujeta a la información disponible en t 1. Normalmente el conjunto de información disponible en el instante t, es decir (t ), es referido a las realizaciones pasadas o rezagos propios del proceso, pero se deja en claro que puede existir una fuente de información adicional relevante. Se de…ne entonces la función de verosimilitud L como el producto de las densidades condicionales parciales, tal que:
T
L =
Yj
f yt t1 , donde 0 es dado.
t=1
Una transformación monótona conocida es la de tomar logaritmos y enseguida maximizar la función Log (L). Como se sabe, dicha transformación no altera el comportamiento esencial de la función objetivo por tanto el punto que maximiza a la transformada es exactamente el mismo que le corresponde a la original. Dicha transformacion permite un trabajo operativo mucho mucho más sencillo, por lo que es frecuentemente utilizada utilizada nuestro ejemplo ejemplo no será la excepción. La funcion de log-verosimilitud es como sigue: T
Log (L) =
X j t=1
Log f yt t1 , donde 0 es dado.
Ahora, es necesario señalar que dicha función de verosimilitud verosimilitud se encuentra asociada al vector de parámetros del proceso generador de datos de yt . Sea el vector de parámetros asociado al proceso generador de datos de y t , el problema se plantea de la siguiente forma: T
max Log (L) =
^
X t=1
Log f yt ; t1 , donde 0 es dado.
j
y por lo tanto = Arg = Arg m ax Log ax Log (L).
64
CAPíTULO 6 APLICACIONES APLICACIONES DE LABORATORIO LABORATORIO
Es importante mencionar que normalmente normalmente las condiciones de primer orden a este problema no se caracterizan por ser lineales, lo que hará necesario el utilizar un algortimo iterativo hasta alcanzar la convergencia necesaria. Es por ello que requerimos un soporte computacional adecuado para llevar a cabo esta tarea, que otro modo no sería posible. A continuación pasaremos revista a los experimentos computacionales que nos van a proporcionar una mejor comprension al tema tratado.
6.3.3. Estimación Estimación por Máxima Máxima Verosim Verosimilitud ilitud en Eviews Eviews El programa Eviews cuenta con una rutina que ejecuta el procedimiento iterativo de encontrar el valor máximo de la función de verosimilitud, dadas la especi…cación de la misma y las condiciones iniciales del ejercicio. En vista de ello, nuestro trabajo se va a enfocar en especi…car de la mejor manera el problema, de modo tal que el resultado sea consistente con los planteamientos iniciales. El objeto logl se encuentra provisto de rutinas generalizas que permiten maximizar una función de verosimilitud respecto a los parámetros parámetros especi…cados del modelo. Lo importante es estar familiarizados con el procedimiento de generar series en el Eviews, puesto que la especi…cación y estimación se basa en un conjunto de series creadas como componentes de la función a maximizar. Supongamos en nuestro caso que hemos identi…cado que una serie sigue un proceso AR(1) y queremos probar empíricamente que esto es cierto. La estimación por máxima verosimilitud nos pide que especi…quemos correctamente la función objetivo creando, en primer lugar, el objeto logl: Logl logl1
La serie que vamos a analizar sigue el proceso:
N (0; (0; 2" )
yt = = y + yt1 + "t , donde " t
Claramente los parámetros a estimar con tres: y , y 2" . Dado el supuesto de normalidad de los errores planteamos la Función de Log-Verosimilitud:
L y ; ; 2" =
T 2
P
log(2 log(2 ) + log 2"
Que en forma resumida resulta:
T
t=1
P
L y ; ; 2" =
T
Lt y ; ; 2"
t=1
(yt y yt 1 )
2 2"
2
6.3 ESTIMACIÓN POR MÁXIMA VEROSIMILITUD
donde: Lt y ; ; 2" = log
(yt y yt 1 )
"
65
1 2 log
2" .
En Eviews vamos a buscar los parámetros iniciales para este ejercicio, siendo los resultados de una regresión MCO los candidatos más adecuados. Estimando la regresión: LS y c AR(1)
Para identi…car mejor los parámetros recomendamos crear vectores de coe…cientes distintos al tradicional. En este caso hacemos: coef(1) mu coef(1) rho coef(1) sigma2
Asignamos los resultados de la regresión a los coe…cientes creados, tomando estos la forma de parámetros iniciales. mu(1)=c(1) rho(1)=c(2) c(3)=mco.@se sigma2=c(3)
Especi…camos la función de verosimilitud tal como fue discutida: logl1.append @logl loglike logl1.append @byeqn logl1.append res=ar1-mu(1)-rho(1)*ar1(-1) logl1.append var=sigma2(1) logl1.append loglike=log(@dnorm(res/@sqrt(var)))-log(var)/2
Eliminamos la primera observación de la muestra, dado que al tomar un rezago esta se pierde: smpl @…rst+1 @last
Con estas especi…caciones el modelo está listo para ser estimado. Por ello, ejecutamos dicha sentencia y observamos los resultados: logl1.ML show logl1
El programa ha hecho lo que le hemos pedido, es decir, una vez que tuvo las especi…caciones sobre la función de verosimilitud no hubo problemas para estimarla. Lo siguiente es evaluar la e…ciencia de este estimador, para lo cual debemos analizar las propiedades estadísticas de los mismos. Un experimento interesante sería generar aleatoriamente muchas series (n) AR(1) que contengan los mismos parámetros, de tal forma que se obtengan los valores MV para cada experimento y se evalúe la convergencia de ellos.
66
CAPíTULO 6 APLICACIONES DE LABORATORIO
6.3.4. Análisis de la e…ciencia del estimador La evaluación de este proceso se hace mediante la obtención de la matriz de covarianzas de los coe…cientes estimados. Luego, es necesario trabajar con las derivadas de esta función de verosimilitud, con el objetivo de construir la matriz de información que se utiliza para comprobar que el estimador sea de mínima varianza. Si la matriz de covarianzas es igual a la inversa de la matriz de información (Cota de Cramer-Rao) entonces se comprueba que el estimador, que asume distribución normal, es de mínima varianza. Capturamos la matriz de covarianzas: matrix cmatr=logl1.@coefcov
Luego construimos la matriz hessiana, aproximada por el producto externo de los gradientes:
logl1.makegrads(n=gradiente) gmu grho gsigma2 matrix g1=@convert(gradiente) matrix g2=@transpose(g1) matrix info=g2*g1 matrix cmatr_info=@inverse(info)
Evaluamos la desigualdad de cramer rao:
6.3 ESTIMACIÓN POR MÁXIMA VEROSIMILITUD
67
matrix cramer_rao=cmatr-cmatr_info
En caso de cumplir la condición, esta matriz debe ser cercana a cero en todos sus elementos. Se recomienda borrar los objetos temporales para un mayor orden en el work…le: d g1 g2 info gradiente gmu grho gsigma2
6.3.5. Test de signi…cancia El test de wald evalúa las distintas restricciones que se pueden imponer a un modelo, evaluando la signi…cancia validez de este conjunto de condiciones. Fijamos la hipótesis nula y evaluamos el resultado de este test: freeze(wald_test) logl1.wald mu(1)=0.5, rho(1)=0.8, sigma2(1)=1 show wald_test
6.3.6. Estimación por Máxima Verosimilitud en Matlab En esta sección estimaremos por Máxima Verosimilitud los parámetros de un proceso ARMA ( p; q ), apuntando a realizar numerosos experimentos computacionales que revelen la distribución de los parametros obtenidos. Simularemos este proceso y luego construiremos la función de verosimilitud para su posterior estimación. Sea el siguiente proceso ARMA (2; 1): yt = 1 yt1 + 2 yt2 + "t 1 "t1 , donde "t N (0; 2" ) Supongamos que 1 = 0;6, 2 = 0;2, 1 = 0;5, " = 1. Luego generamos un proceso de 1000 observaciones: % Generamos un proceso ARMA(p,q) c=0;% Intercepto A=[1 -0.6 0.2]; % Componente AR B=[1 0.5];% Componente MA T=1000;% Número de Obs. sigma=1;% Desv. Estándar del término de perturbación p=length(A); q=length(B); [y e]=armapq(A,B,T+p,sigma); y=c+y;
Antes de seguir con el procedimiento es necesario crear una función de verosimilitud a maximizar. De preferencia se genera un archivo .m adicional para luego llamar a una función de optimización que trabaje con el vector de parametros objetivo.
68
CAPíTULO 6 APLICACIONES DE LABORATORIO
Creamos la función en base al proceso ARMA (2; 1) generado. Se asume distribucion normal para el término de perturbación:
P T
L (1 ; 2 ; 1 ; " ) =
Log
t 1
yt 1 yt
1
y
2 t2
"
son dados.
+1 "t
1
, donde y 0 ; y1 ; "0
function L=LogL(y,e,theta,p) % Se requiere de preferencia p >q % p: Términos AR n=length(theta); T=length(y); % T+p en este caso AR=theta(1:p); for i=1:p ARi(i)=AR(p-i+1); end MA=theta(p+1:n-1); q=length(MA); for i=1:q MAi(i)=MA(q-i+1); end sigma=theta(n); L=0; for i=p+1:T lagsy=y(i-p:i-1); lagse=e(i-p:i-p+q-1); L=L+log(normpdf(y(i,1),ARi*lagsy+MAi*lagse,sigma)); end L=-L;
Como podemos apreciar, estamos tomando el negativo de la funcion de verosimilitud. Esto es desde que la herramienta de Matlab permite minimizar funciones de varias variables. En este caso, al tomar el negativo, estaríamos maximizando la misma. Creamos una función adicional para que solamente quede como argumento: function Lf=LogL_F(theta) global y e p Lf=LogL(y,e,theta,p);
Con la función de verosimilitud correctamente construida pasamos a realizar el experimento:
6.3 ESTIMACIÓN POR MÁXIMA VEROSIMILITUD
69
% Estimación por Máxima Verosimilitud de los parámetros de un ARMA(p,q) % % Fernando Pérez Forero % 11/12/05 clear % Generamos un proceso ARMA(p,q) c=0;% Intercepto A=[1 -0.3 0.2]; % Componente AR B=[1 0.3];% Componente MA T=1000;% Número de Obs. sigma=1;% Desv. Estándar del término de perturbación p=length(A); q=length(B); [y e]=armapq(A,B,T,sigma); y=c+y; AR=A(2:p); AR=-1*AR; p=length(AR); MA=B(2:q); MA=-1*MA; q=length(MA); y0=randn(p,1); y=[y0;y]; e0=randn(q,1); e=[e0;e]; AR_CI=[0.29 -0.19]; MA_CI=[-0.29]; sigma_CI=0.99; theta=[AR_CI,MA_CI,sigma_CI]; global y e p % Seleccionando el vector de opciones default options = optimset(’fminsearch’); % Ajustando las opciones options.MaxFunEvals= 30000; options.LargeScale=’o¤’; options.MaxIter= 15000; [theta_…n,LogF,exit‡ag,output] =fminsearch(’LogL_F’,theta,options); LogF=-LogF; theta_…n
70
CAPíTULO 6 APLICACIONES DE LABORATORIO exit‡ag
^
Como resultado obtenemos el vector que maximiza la función de logverosimilitud:
theta_…n = [0.5726, -0.3259, -0.0181, 0.9512] La agenda pendiente se centra en la obtención de las varianzas de estos estimadores, así como el realizar experimentos múltiples para veri…car la consistencia de los mismos.
6.3.7. Asignación Primera parte En la seeción II se estimaron los parámetros especi…cando una función de verosimilitud. Ahora, se sabe bien que en el caso de los AR(1) la primera observación tiene una regla de correspondencia distinta, lo cual hace que la especi…cación y estimación sea condicional a qué observación se está utilizando. La regla de correspondencia completa sería: f (y; ) =
8>< p ( ) >: 1 2(1 2 )
1 exp 2
p
exp
yt
2
2
y 1
2
, t = 1
12
2
(yt y yt 1 )
2( 2 )
, t > 1
En base al código anterior, modi…que la especi…cación para satisfacer esta restricción. Pistas: Utilice una variable dummy que distinga los dos intervalos de la regla de correspondencia. La función recode ejecuta sentencias en base a una condición inicial, de no cumplirse ejecuta otra sentencia indicada.
Segunda parte Evalúe la consistencia de las muestras grandes generando 500 experimentos (!Repl) tanto para un ejercicio de 35 como de 1000 observaciones, utilizando de preferencia el resultado de la pregunta anterior. Almacene los resultados de cada replicación en un vector columna, de forma que se pueda evaluar la distribución asintótica de cada estimador. Si es posible, haga lo mismo con las desviaciones estándar. Pistas: Para generar un AR(1) ejecute los siguientes comandos: work…le u 1 !Obs
6.3 ESTIMACIÓN POR MÁXIMA VEROSIMILITUD
71
!Obs=500 !mu_pobl=0.5 !rho_pobl=0.8 series AR1=0 smpl @…rst+1 @last series AR1=!mu_pobl +!rho_pobl *AR1(-1)+nrnd smpl @all
Para crear los vectores de resultados ejecute los siguientes comandos: vector(!Repl) mudist, Vector de parámetros mu estimados. vector(!Repl) stdmu, Vector de desv. de mu de cada experimento. vector(!Repl) rhodist , Vector de parámetros rho estimados. vector(!Repl) stdrho, Vector de desv. de rho de cada experimento. vector(!Repl) sigma2dist, Vector de parámetros sigma2 estimados. vector(!Repl) stdsigma2 , Vector de desv. de sigma2 de cada experi-
mento.
Para almacenar el resultado de cada parámetro del experimento i en el
vector ejecutar:
mudist(!i)=mu(1) , Almacenamos los resultados de la simulación i para
mu.
72
CAPíTULO 6 APLICACIONES DE LABORATORIO
6.4. Ejercicios 1. Bajo que restricciones sobre los parámetros del siguiente modelo, este es estacionario en media y varianza yt = 1 yt1 + 2 yt2 + t
donde , t N (0; 2 ): 2. Asumiendo que los parámetros del modelo anterior cumplen con las restricciones para garantizar estacionariedad, calcular el la respuesta de yt h periodos hacia adelante ante un choque en t , esto es calcular @y t+h = @ t
Suponga, que este modelos se puede expresar como (1
1L) (1 2L) yt = t
Escriba su respuesta en función de 1 y 2 :
3. Explique cuales son las propiedades del estimador de máxima verosimilitud 4. Calcule la predicción el error de predicción y intervalos de con…anza para el siguiente modelo AR(2) yt = 1 yt1 + 2 yt2 + t
5. Demuestre que la esperanza condicional en t es el predictor lineal óptimo. 6. Es el siguiente proceso M A(2) estacionario en covarianzas? Y t = (1 + 2;4L + 0;8L2 ) t
donde, E ( t t+s ) =
8< :
1 para s = 0 0 en otro caso
si este es el caso, calcule sus autocovarianzas
73
6.4 EJERCICIOS 7. Es el siguiente proceso AR(2) estacionario? (1
donde, E ( t t+s ) =
8< :
1;1L + 0;18L2)Y t = t 1 para s = 0
0 en otro caso
si lo es, calcule sus autocovarianzas. 8. Dado el siguiente modelo AR(2) estacionario en covarianzas yt = 1 yt1 + 2 yt2 + t
donde , t N (0; 2 ): Encuentre su representación en medias móviles. ( Hint: use la siguiente condición: B(L) A(L) = 1
donde:
B(L) = 1 1 L 2 L2 A(L) = 1 + 1 L + 2 L2 + 3 L3 ::::::
9. Dado el siguiente modelo AR(P), yt = 1 yt1 + 2 yt2 :::: p yt p + t
donde , t N (0; 2 ) Escríbalo en forma de estado espacio, y especi…que bajo qué condiciones es estacionario. : 10. Dado el siguiente modelo AR(1) yt = yt1 + t
donde , t N (0; 2 ) Escriba la seudo función de verosimilitud, encuentre el estimador de MV de y de 2 , asi como de sus respectivas varianzas. 11. Muestre que la precisión de la proyección de un modelo AR(1) disminuye a medida que el horizonte de predicción aumenta.
74
CAPíTULO 6 APLICACIONES DE LABORATORIO
12. Explique en que consisten las siguientes pruebas estadísicas: a ) La prueba de normalidad de Jarque-Bera b ) Prueba de autocorrelación de LM c ) Prueba de heterocedasticidad de LM d ) Prueba de heterocesdacidad de White, e ) La prueba de razón de verosimilitud.
13. Explique en detalle cuál es el procedimiento para identi…car apropiadamente un modelo ARIMA(p,d,q). 14. Explique cuales de las siguientes proposiciones son verdaderas o falsas a ) Todo proceso estacionario ARMA(p,q) tiene una representación de medias móviles M A ( )
1
b ) La matriz de varianza covarianzas del estimador de máxima verosimil-
itud esta dado por el producto orthogonal de la gradiente de la función de verosimilitud respecto a los parámetros. c ) En el modelo AR(2) siguiente yt = 1 yt1 + 2 yt2 + t
la respuesta de largo plazo de yt ante un choque en t es 11 1
2
Capítulo 7
Vectores Autorregresivos1 Los vectores autorregresivos son la generalización de los modelos autorregresivos univariantes al caso multivariante, tienen la ventaja de incorporar mucha mas información, esto es para la proyección de cada variable endogena, se considera toda su historia, y también la historia de todas las otras variables que forman parte del modelo VAR(p). Asimismo, se consideran varias fuentes de choques como las generadoras de ‡uctuaciones.Veamos el siguiente ejemplo:
24 35 24 y t
t
=
11 12 21 22
35 24 35 24 35 y t1
+
t1
v1;t v2;t
Donde yt y t representa la brecha producto y la in‡ación. Hay varias características de este modelo que vale la pena resaltar, primero, la dinámica tanto de la in‡ación como del producto dependen ahora de los valores pasados tanto de la in‡ación como de la brecha producto, esto es, el modelo para proyectar yt esta dado por. yt = 11 yt1 + 12 t1 + v1;t
Segundo, cada variable depende de todas las fuentes de ‡uctuaciones puesto que v1;t es una combinación lineal de las fuentes de choques primarias en el modelo. Esta información se obtiene del hecho que 1
La mayoría del material discutido en estas notas es un resumen de los capítulos 10 y 11 de Hamilton (1994) y del libro de Amisano y Giannini "Topics in Structural VAR Econometrics"
75
CAPíTULO 7 VECTORES AUTORREGRESIVOS2
76
E (V t V t0 ) =
24
11 12 12 22
35
La covarianza entre el choque 1 y 2 no es cero, es decir estas variables no son independientes. Esto quiere decir que si v 2;t cambia, también lo hará v1;t afectando la brecha producto. Es por ello que antes de realizar ejercicios de simulación dinámica se deben ortogonalizar estos residuos. Cuando los residuos no están ortogonalizados la representación del modelo VAR se denomina forma reducida. Veamos a continuación un ejemplo como un modelo dinámico estructural se puede escribir como un VAR(1). Un modelo típico de una economía cerrada sin expectativas yt = 1 yt1 2 (it (t )) + gt t = t1 + y t + t it = ' t + 'y yt + t
donde fyt ; t ; it g representan las variables endógenas del modelo, y fgt ; t ; it g las variables exógenas. La primera ecuación representa una curva IS dinámica, la segunda una curva de Phillips, y la tercera una regla de Taylor. Escribiendo este modelo en forma matricial tenemos:
26 64
1
a2
a2
1
0
'y
'
1
37 26 75 64
yt t it
37 75
=
26 64
1 0 0 0 0 0 0 0
37 26 75 64
yt1 t1 it1
37 75
+
26 64
gt t t
37 75
de forma compacta, el sistema anterior puede escribirse como: B0 Z t = B1 Z t1 + t
donde Z t =
B1 =
26 64 26 64
yt
37 75
1
t ; B0 = it
1 0 0
26 64
37 75
a2
a2
1
0
'y
' 1
0 0 ; t = 0 0 0
26 37 64 75 gt
37 75
;
t t
Multiplicando ambos lados de esta ecuación por B 01 tenemos que:
7.1 ESTACIONARIEDAD Y MOMENTOS
77
Z t = B01 B1 Z t1 + B01 t
y haciendo 1 = B01 B1 y V t = B01 t notamos que el modelo estructural dinámico anterior puede escribirse como un modelo VAR(1). El modelo es estructural porque impone restricciones sobre el numero de variables que afectan cada ecuación y sobre la relación que existe entre los diversos choques. En el caso particular de este ejemplo, el choque de gasto publico, de oferta y el choque monetario se asume que son ortogonales entre si, razón por lo cual su matriz de varianza covarianza es escalar. Sin embargo, la forma reducida del modelo VAR(1), la que se obtiene al multiplicar el sistema por B01 tiene choques que no son ortogonales, V t , su matriz de varianza covarianza no es escalar. E (V t V t0 ) = B01 E ( t 0t )B010
Sin embargo, las correlaciones que se observan en los choques de la representación en forma reducida contienen información sobre las relaciones contemporáneas de las variables bajo análisis, en este caso in‡ación y producto. Para poder realizar simulaciones dinámicas con este modelo, primero debemos estimar B 0 , puesto que dada esta información podemos trasformar los residuos de la forma reducida, que son estimables directamente, en residuos ortogonales que pueden utilizarse para hacer simulaciones dinámicas. t = B0 V t
En lo que sigue vamos a revisar algunos conceptos básicos sobre VAR, que son la generalización de lo que hemos visto para modelos univariantes, esto es, el calculo de los primeros y segundos momentos de un VAR(p), condiciones de estacionariedad, el procedimiento para invertir in VAR(p), el procedimiento de estimación, pruebas de hipótesis, y predicción, …nalmente regresaremos al tema de simulaciones dinámicas para analizar las funciones de impulso respuesta, la descomposición de varianza, y varias formas de estimar t , es decir VAR(p) estructurales.
7.1. Estacionariedad y momentos Un VAR(p) típico tiene la siguiente forma: Y t = o + 1 Y t1 + 2 Y t2 :::: + p Y t p + t
CAPíTULO 7 VECTORES AUTORREGRESIVOS3
78
donde Y t y t son vectores de dimensión k 1, y j matrices de orden k k que contienen los parámetros del modelo, y la matriz de varianza covarianza de t satisface: para t = s 0 E ( t t+s ) =
8< :
0 en otro caso
Al igual que en el caso de modelos univariantes, para garantizar las estacionariedad débil de estos modelos se tiene que restringir sus primeros y segundos momentos a ser independientes del tiempo. Bajo el supuesto de estacionariedad. E (Y t ) = (I
1 2 :: p)1 o
Asimismo, utilizando el operador de rezagos, el modelo anterior puede escribirse como: (L)Y t = t
donde, el polinomio de rezagos se de…ne como: (L) = I
1L 2L2 :: pL p
El modelo anterior puede escribirse también en forma estado espacio, es decir como un VAR(1) de la siguiente forma: Z t = F Z t1 + V t
Donde estas nuevas variables se de…nen a continuación como:
Z t(kp1) =
F (kpkp) =
y
26 66 4 26 66 664
Y t
o Y t1 o .. .
Y t p1 1 I k 0
.. . 0
37 77 5
o 2 2 0 0 I k 0 .. .. . . 0 0
p1 p 0
0
0
0
.. .
.. .
I k
0
37 77 775
79
7.1 ESTACIONARIEDAD Y MOMENTOS
V t(kp1) =
26 66 4
t 0
.. . 0
37 77 5
donde la matriz de varianza covarianza de este nuevo vector de errores esta dado por: Q para t = s 0 E ( t t+s ) =
donde
Qkpkp =
26 66 664
8< :
0 0
.. . 0
0 en otro caso
0 0 0 0 0 0
.. .. . . 0 0
37 77 775
Las condiciones de estacionariedad del modelo AR(p) pueden encontrarse de manera sencilla, a partir de la matriz F. Se dice un modelo VAR(p) es estacionario si sus primeros y segundos momentos son independientes del tiempo. Esta condición se satisface si la matriz F , correspondiente a la representación estado espacio del VAR tiene todos los valores propios i tales que ji j < 1: Los segundos momentos del VAR(p), entonces pueden calcularse a partir de:
26 66 4 0
E Z t Z t = =
0 0
1
.. .
1
0 p2
p1 0
p2 0
37 77 5
.. . donde, j representa la autocovarianza de orden j de Y t . Para calcular, procedemos como el caso de un ar(1).
0
p1
0
0 = E (F Z t1 + V t ) (F Z t1 + V t )0 = F 0 F 0 + Q
De igual manera, podemos calcular:
(7.1)
80
CAPíTULO 7 VECTORES AUTORREGRESIVOS4
E Z t Z t1 = E (F Z t1 + V t ) (Z t1 )0 = 1 0
este procedimiento puede generalizarse para obtener: 0
E Z t Z t j = F j 0 = F j 1
Notese, que las formulas que estamos derivando, son exactamente las mismas que las derivadas en clases anteriores para el caso de un AR(1), si asumimos que F ,Q y 0 son escalares. Para obtener una forma cerrada de la solución de 0 debemos hacer uso de dos operadores adicionales, el operador vec y el producto Kronecker. El operador vec convierte, una matriz de n por m, en un vector nm por 1, apila las columnas de la matriz una debajo de otra en el vector. Por ejemplo. dada la matriz, A =
el vec(A) =
26 66 4
a11 a21 a12 a22
37 77 5
24
a11 a12 a21 a22
35
Proposition 2 El operador vec satisface la siguiente propiedad: vec(ABC ) = (C 0 A) vec(B) donde el símbolo denota el producto kronecker. Aplicando el operador vec a la ecuación 7.1 tenemos que: vec(0 ) = vec(F 0 F 0 ) + vec(Q) vec(0 ) = F
Fvec(0) + vec(Q)
vec(0 ) = [I r
F F ]1 vec(Q)
donde r = (kp)2 si la matriz F tiene todos sus valores propios estacionarios, este también sera el caso de la matriz F F entonces se puede garantizar que existe 0 .
7.2 CAUSALIDAD EN EL SENTIDO DE GRANGER
81
7.2. Causalidad en el sentido de Granger Una pregunta relevante en el análisis VAR es el referido al contenido de información de una variable respecto a otra. Para responder esa pregunta utilizamos el concepto de causalidad en el sentido de Granger, que no se re…ere a una relación causa efecto, sino mas bien es un concepto de precedencia estadística. Se dice que una variable xt no causa otra variable yt en el sentido de Granger, si xt no ayuda a predecir yt , mas formalmente, si EC M (yt+s =I t ; xt ) = ECM (yt+s =I t ) donde I t = fyt1: yt2 ::g Para probar estadísticamente la causalidad en el sentido de Granger, se puede utilizar pruebas estadísticas F sobre los coe…cientes de un modelo VAR(p). Por ejemplo en el modelo VAR de…nido anteriormente. para la in‡ación y el producto:
24 35 24 35 24 35 24 35 24 35 24 35 y t
=
t
p11 p12 p21 p22
111 112
y t1
121
t1
122
y t p
t p
:: +
+ ::
v1;t v2;t
La in‡ación no causa en el sentido de Granger al producto si, J 12 = 0, para todo j:Esto es la in‡ación es no informativa de la brecha producto.
7.3. La representación en medias móviles Todo VAR(p) estacionario tiene una representación equivalente como un V MA(1), la representación de medias móviles es muy conveniente puesto que permite fácilmente calcular la predicción en un VAR(p), los impulsos respuesta y la descomposición de varianza. Dado el siguiente modelo VAR(p) estacionario Y t = o + 1 Y t1 + 2 Y t2 :::: + p Y t p + t
Existe su representación en medias móviles dada por:
Y t = 0 + t + 1 t1 + 2 t2 + :::::
Para invertir el modelo AR(p), tenemos que encontrar las matrices j como función de las matrices j . Para ello, expresamos ambos modelos utilizando el operador de rezagos: (L)Y t = t
CAPíTULO 7 VECTORES AUTORREGRESIVOS5
82
Y t = (L) t
donde debe cumplirse que: (L)1 = (L)
entonces: (L)(L) = I
I
1 L
2L2:::: pL p
I + 1 L + 2 L2 + ::: = I
Desarrollando esta ecuación e identi…cando términos tenemos las siguiente secuencia recursiva para los j como función de los j 1 = 1
2 = 1 1 + 2
s = 1 s1 + 2 s2 ::: p s p
Notese, que una vez obtenida la version de medias móviles del VAR(p) en forma reducida, se puede obtener la version de medias móviles en la forma structural, notando que: t = B01 V t
Entonces, el VAR( 1) en forma estructural esta dado por: Y t = (L)B01 V t
Dada la representación en medias móviles es muy sencillo obtener la función de impulso respuesta y la descomposición de varianza del modelo.
7.4 FUNCIÓN IMPULSO RESPUESTA
83
7.4. Función impulso respuesta Dada la representación de medias móviles de la forma reducida Y t = c + t + 1 t1 + 2 t2 + :::
La representación de medias móviles de la forma estructural esta dada por: Y t = c + B01 V t + 1 B01 V t1 + 2 B01 V t1 + :::
Entonces, la función impulso respuesta describe la respuesta de todo el vector de variables Y t+s a un impulso de una sola vez en y j;t , manteniendo todas las demás variables constantes. Analíticamente, el impulso respuesta se mide como: @Y t+s = s b j @v j:t
Donde:
h i
B01 = b1 b2 : : : bk
En el ejemplo del modelo estructural discutido anteriormente, se pueden calcular los impulsos respuesta de la brecha producto, la in‡ación y las tasas de interés a un choque de gasto publico, costos, y de política monetaria. En cada caso, la función impulso respuesta muestra la respuesta temporal de estas tres variables endogenas a los choques:
7.5. Descomposición de Varianza Puesto que existe multiples fuentes de ‡uctuaciones en un modelo VAR(p), cada una correspondiente a un choque estructural, es importante medir que tan importante es cada fuente de ‡uctuaciones explicando la volatilidad total de cada serie endogena. Por ejemplo, nos gustaría saber de los 3 choques estructurales cual es el mas importante para explicar el producto o la in‡ación, son los choques de demanda o los de oferta. Este es precisamente la información que nos brinda la descomposición de varianza. Nuevamente, utilizando la representación de medias móviles del modelo estructural Y t = c + B01 V t + 1 B01 V t1 + 2 B01 V t1 + :::
CAPíTULO 7 VECTORES AUTORREGRESIVOS6
84
El error de predicción, s periodos hacia adelante esta dado por: Y t+s
Y t+s=t = B01V t+s + 1B01V t+s1 + 2B01V t1 + :::s1B01V t+1
La varianza del error de predicción esta dada por:
E (Y t+s
Y t+s=t)(Y t+s Y t+s=t)0
donde
= + 1 01 + :::s1 0s1
E (V t V t0 ) = E ( t 0t ) = = B 01 B010
Puesto que t = B 01 V t = b 1 v1;t + b2 v2;t : + :bk vk;t = b 1 b10 var(v1;t ) + b2 b20 var(v2;t )::: + bk bk0 var(vk;t )
La varianza del error de predicción puede escribirse como la suma de la contribución marginal a la varianza total de la varianza de cada choque: j=k
X
E (Y t+s Y t+s=t )(Y t+s Y t+s=t )0 =
j=1
var(v j;t ) b j b j 0 + 1 b j b j 0 01 + :::s1 b j b j 0 0s1
Entonces, la contribución marginal del choque j en la varianza del error de predicción de la variable i esta dado por: 100
var(v j;t)diag
b j b j 0 + 1 b j b j 0 01 + :::s1 b j b j 0 0s1 (i; i)
A medida que el horizonte de predicción se hace mayor, la varianza del error de predicción converge a la varianza no condicional de la serie 0 de manera análoga a lo que demostramos en la clase de predicción para el caso de modelos univariantes.
85
7.6 ESTIMACIÓN DE MODELOS VAR(P)
7.6. Estimación de modelos VAR(p) La estimación de un modelo VAR(p) se hace utilizando la seudo función de verosimilitud, esto es es utilizando solamente las distribuciones de probabilidad condicionales. Como se discutió en clases anteriores en este caso los estimadores de maxima verosimilitud coinciden con los de mínimos cuadrados ordinarios.Dado el siguiente modelo VAR(p) estacionario Y t = o + 1 Y t1 + 2 Y t2 :::: + p Y t p + t
con t N (0; ) tenemos que: Y t Y t1 Y t2 :::Y t p+1 :
j
j
N (o + 1Y t1 + 2Y t2:::: + pY t p; )
El modelo anterior se puede escribir de forma compacta apilando las variables exogenas, en este caso predeterminadas en un vector de variables X t tal que:
X t =
26 66 4
1 Y t1
.. .
Y t p
37 77 5
y = [o 1 ::: p ]
entonces, el VAR(p) puede escribirse como: Y t = 0 X t + t
Por tanto la función de verosimilitud estaría dada por:
L() =
X h
Tk T log(2) + log 1 2 2
T
1 2 t=1
i
0 Y t 0 X t 1 Y t 0 X t
Esta función de verosimilitud coincide con la de un modelo lineal general, por tanto el estimador de maxima verosimilitud de 0 y esta dado por:
"X #"X # T
0k(kp+1) =
Asimismo se tiene que:
T
Y t X t0
t=1
X t X t0
t=1 t
1
CAPíTULO 7 VECTORES AUTORREGRESIVOS7
86
"X #"X # T
j10(kp+1) =
T
y j;t X t0
X t X t0
t=1
1
t=1 t
Este estimador equivale a estimar por mínimos cuadrados ordinarios y j;t en X t0 . De manera similar,
" # X b T
1 = T
0
t t
t=1
Bajo el supuesto que t es independiente e idénticamente distribuido con media cero, y varianza se cumple que ( proposición 11.1 Hamilton (1994): 1.
1 T
P T
t=1t
X t X t0
1
p ! Q
p
b bp ! b ! !
2. ! 3. 4.
p
T T
d
N (0;
Q)
Es decir, bajo el supuesto de independencia, idéntica distribución y homocedasticidad y no correlación serial, el estimador de maxima verosimilitud, que coincide con el estimador de MCO es consistente y tiene una distribución asintótica que es normal. Esto resultado permite hacer pruebas de hipótesis sobre los parámetros del modelo de manera muy sencilla. Así por ejemplo podemos imponer un conjunto de restricciones lineales sobre de la siguiente forma: R = r
Que dado el supuesto de normalidad de implica que:
p
b ! !
T RT
r
d
N (0; R (
Q) R0)
Lo que nos permite de…nir el siguiente estadístico de Wald,
b ! 0@ 0@ "X # 1A 1A b ! T
W = RT
r
R
X t X t0
t=1 t
1
1
R0
RT
r
7.7 IDENTIFICACIÓN DE MODELOS VAR(P)
87
0 Donde W 2 (m) donde m es el numero de restricciones que se imponen. De igual manera, es sencillo implementar pruebas de hipótesis utilizando el estadístico de la razón de verosimilitud. RV = 2(L1
L0) = T
1 log 0
1 log 1
Este estadístico tiene una distribución 2 (m) donde m es el numero de restricciones que se imponen
7.7. Identi…cación de modelos VAR(p) De manera similar al caso de modelos univariantes con modelos multivariantes de series de tiempo se requiere también la correcta identi…cación del orden de VAR. Para ello seguimos el mismo principio que en el caso de modelos univariantes, esto es el p es aquel que genera un vector de residuos t bien comportado, esto es sin presencia de correlación serial, heterocesdacidad, y que se comporte siguiente una distribución normal. Por tanto para determinar p se sigue un procedimiento secuencial, se estima el VAR con un p inicial, p0 , luego se evalúa la calidad del residuo, si este no presenta las características deseadas se aumenta el p hasta encontrar el valor adecuado. Para evaluar los residuos, además de pruebas multivariantes de correlación serial, heterocesdacidad y normalidad, se utilizan criterios de información multivariantes, y pruebas de hipótesis para determinar la signi…cancia de rezagos adicionales. Los criterios de información mas conocidos son el de Akaike, AIC, el bayesiano-Schwarz-Bayesian y Hannan-Quinn (HQ):
AIC ( p) = log 1 + BI C ( p) = ln 1 + HQ( p) = ln 1 +
2 pk2 T
ln(T ) pk2 T
2 ln(ln(T )) pk2 T
El criterio de AIC sobreestima el verdadero valor de p asintóticamente, mientras que los BIC y HQ son asintóticamente consistentes, si es que el p < p m ax. Para realizar la prueba de normalidad de los residuos se puede utilizar un extension al caso multivariante de la prueba de Jarque Bera, esta prueba
CAPíTULO 7 VECTORES AUTORREGRESIVOS8
88
compara el tercer y cuarto momento de la distribución con los valores teóricos de estos momentos para el caso de una normal. Esta prueba se aplica a los residuos ortogonalizados. La ortogonalización de los residuos se puede alcanzar por cualquiera de los métodos que vamos a discutir en la sección siguiente. Dado un vector de residuos de la forma estructural tal que: V t = B 0 t
N (0; I k )
Entonces, tenemos que:
2 3 2 p 4 5 4 T
m3
m4
N (0;
3
6I k
0
0 6I k
35
)
donde m3 y m4 representa el momento de tercer y cuarto orden. Para las pruebas individuales de normalidad, se utiliza el estadístico típico de Jarque Bera, m23 (m4 3)2 J = ( + ) 6 24
2(2)
Para la prueba de normalidad conjunta se realizan las siguientes pruebas. 3 = 4 =
T (m4
T m03 m3 6
2(k)
3)0 (m3 3) 2(k) 24
2(2k)
= 3 + 4
Capítulo 8
VAR Estructurales En el capítulo anterior, para derivar las ecuaciones correspondientes al impulso respuesta y a la descomposición de varianza, se asumió que la relación entre errores de la forma reducida y de la forma estructural era conocida, esto es, se suposo que se conocía la matriz, B 0 , y por lo tanto era posible determinar, a partir de los residuos de la forma reducida, iguiente relación
V t = B0 t
Donde E t t 0t = y E t (V t V t0 ) = Sin embargo, en la practica se estima de la representación de forma reducida del VAR, y en general las matrices B0 y no están identi…cadas, esto es no es posible determinar de manera única estas dos matrices utilizando la información contenida en la estimación de la forma reducida. Puesto que es una matriz simétrica semide…nida positiva, solo contiene k(k+1) condiciones independientes que se pueden utilizar para estimar 2 los k2parámetros de B0, se requieren entonces con k(k21) restricciones adicionales para poder identi…car los coe…cientes de la matriz de relaciones contemporáneas de las variables. Estas restricciones adicionales se pueden imponer de varias maneras, en lo que sigue vamos a analizar tres formas alternativas de alcanzar la identi…cación de los parámetros de la forma estructural. Cuando el modelo estructural esta identi…cado se puede obtener de manera sencilla los parámetros de la forma estructural a partir de los parámetros de la forma reducida. Para que un conjunto de parámetros estructurales este identi…cado debe satisfacer dos tipos de condiciones simultáneamente, las condiciones de orden y de rango. La condición de orden se re…ere al numero 89
90
CAPíTULO 8 VAR ESTRUCTURALES
de restricciones adicionales que se deben imponer en el modelo para alcanzar su identi…cación, esto es se requiere por lo menos k(k21) restricciones adicionales. Sin embargo, aun si se cumple con esta condición, es posible que los parámetros de la forma estructural no puedan identi…carse, se requiere además que se satisfaga la condición de rango. Para entender las implicancias de la condición de rango, veamos como se obtienen los parámetros estructurales del modelo VAR utilizando el estimador de maxima verosimilitud. La función de verosimilitud del VAR estructural esta dada por:
i
T2k log(2) T 2 log B01B01 T 1 (Y 0 X )0 1 (Y t 0 X t ) B01 B01 t t 2
L =
Ph
t=1
Dada la estimación de , la función de verosimilitud anterior, puede expresarse de manera compacta como :
L =
T2k log(2) + T 2 log jB0j2 T 2 log jj2 T 2 traza
h
B01 B01
bi 1
Uno puede demostrar, que la función de verosimilitud anterior alcanza su máximo respecto a las matrices , B0, si existen matrices únicas que satisfacen el siguiente sistema de ecuaciones no lineales:
b
B01 B01 =
Utilizando el operador vech al sistema de ecuaciones, tenemos que:
b
V ech() = vech
B01 B01
Para que exista una solución única de este sistema se requiere que matriz hessiana de V ech() respecto a los parámetros que se buscan, tengan el numero de columnas linealmente independientes igual al numero de parámetros que se buscan estimar. Para veri…car esta condición en la practica, el procedimiento mas sencillo es calcular numéricamente el rango de la matriz hessiana de V ech() utilizando un supuesto sobre los parámetros a estimar, y compararlo con el numero de parámetros que se están estimando. Si la condición se satisface se ha encontrado una potencial solución de lo contrario se tienen que modi…car las restricciones de identi…cación apropiadamente hasta lograr la identi…cación de los parámetros .
b
b
91
8.1 DESCOMPOSICIÓN DE CHOLESKI
8.1. 8.1. Desco Descompo mposic sición ión de Chole Choleski ski La forma mas sencilla de alcanzar identi…cación es utilizar la descomposición de choleski, que impone una estructura recursiva en el modelo. Toda matriz simétrica de…nida positiva puede descomponerse de la siguiente forma: = P P 0
Donde P es una matriz triangular inferior. Entonces dado P , y el vector de residuos de la forma reducida, el vector de residuos estructurales esta dado por: V t = P = P 1 t
Lo que implica que
B0 = P = P 1
Supongamos del siguiente ejemplo en un VAR(2), en este caso la estimación de la matriz de varianza covarianza de los residuos nos dan 3 condiciones para identi…car 4 parámetros, los correspondientes a B 0 necesitamos por lo menos una condición adicional para poder identi…car el modelo. =
24
11 12 21 22
35
24 35 24 35 24 35 24 35 24 p p11 0
En este caso, se busca la matriz P =
p21 p22
De la condición anterior p11 0
p21 p22
de donde
p11 p21 0 p22
=
2 p11
p11 p21
2 + p 2 p11 p21 p21 + p22
p11 =
p21 =
p22 =
11
p 21
11
s p 22
21 11
=
11 12 21 22
35
92
CAPíTULO 8 VAR ESTRUCTURALES ESTRUCTURALES
El sistema se resuelve de manera recursiva. recursiva. Notese que la descomposición de Choleski impone una estructura recursiva, esto es, la primera variable en el sistema es la mas exogena, es contemporaneamente independiente de las otras variables, esta determinada únicamente por su propio choque. La segunda depende de la primera y su propio choque, la tercera de la primera y la segunda y así sucesivamente.
24 35 24 35 24 35 v1;t
=
v2;t
p11 0
1;t
p21 p22
2;t
Esta estructura recursiva puede que no sea la mas apropiada desde el punto de vista teórico, por lo que en muchos casos se requieren estrategias de identi…cación identi…cación alternativa alternativas. s.
8.2. Restricci Restriccione oness contempo contemporáne ráneas as Una forma alternativa de alcanzar la identi…cación del modelo es imponer directamente en la matriz B 0 y restricciones de exclusion, esto es hacer que algunas de los elementos sean cero o satisfagan restricciones lineales. Este tipo de restricciones de identi…cación han sido utilizadas para identi…car choques de política monetaria en economías cerradas y abiertas en varios trabajos de investigación. Por ejemplo Bernanke y Mihov, Christiano Eichenbaum y Evans, Leeper, Sims y Zha.
8.2.1. 8.2.1. Bernan Bernanke ke y Miho Mihov A continuación presentamos dos ejemplos de identi…cación de vectores autoregresivos autoregresivos estructurales utilizando restricciones sobre los parámetros de la matriz contemporanea Bo . El primero corresponde al trabajo de Bernanke y Mihov (1998), en el que se utiliza un modelo del mercado interbancario para identi…car identi…car choques choques de política política monetaria monetaria.. El trabajo utiliza utiliza datos datos trimestrales para los EE.UUde 1965 a 1996. El modelo VAR que estiman comprende dos grupos de ecuaciones, las primeras que contienen la dinámica de las variables que no son de política, Y t y las variables de política, P t . El modelo puede escribirse entonces como: k
Y t =
X i=0
Bi Y ti +
k
X i=0
C i P ti + Ay V t
8.2 RESTRICCIONES CONTEMPORÁNEAS k
P t =
k
X i=0
93
Di Y ti +
X i=0
Gi P ti + AP V tP
El conjunto de variables Y t esta dado por, el producto bruto interno, el de‡actor implícito del GDP, y un índice de precios de comodities Y t= [GDP t ; t ; W P It ]
El conjunto de variables de política esta conformado por la tasa de interés de política, F F R, el nivel de reservas no prestadas, N BR y el nivel de reservas totales j.T R P t = [FFR;NBR;TR] FFR;NBR;TR]
El nivel de reservas prestadas corresponde a las reservas generadas a través de la ventanilla de descuento del banco central, mientras que las no prestadas, a las generadas por operaciones de mercado abierto realizadas por el banco central. Para indenti…car los parámetros de la forma estructural se tienen que estimar k2 = 36 parámetros de la matriz B01 . k(k2+1) = 21 condiciones se obtienen de la matriz de varianza covarianza de los residuos
b
= B01 B01
Por lo tanto se requieren de 15 restricciones adicionales para identi…car exactamente el modelo. Para ello se imponen el siguiente grupo de restricciones: 1. Las varia variables bles de política política no tienen efectos efectos contempor contemporáneo áneoss sobre el resto del sistema, C 0 = 0: Puesto que el sistema incluye seis variables, tres de política y tres de no política, la restricción anterior impone 9 restricciones adicionales de identi…cación 2. El subsistema subsistema de variabl variables es de política política satisface satisface las restricci restricciones ones impuestas por un modelo simple del mercado de dinero. Este impone 5 restricciones restricciones adicionales. adicionales. Considerando, los residuos de la forma reducida y forma estructural del subsistema de variables de política monetaria, el modelo puede escribirse como: "T R =
"F F R + vd
94
CAPíTULO 8 VAR ESTRUCTURALES "BR = " T R
"NB R = ("F F R "DI S ) + vB
"N BR = d vd + B v B + v S
La matriz B01 se determina resolviendo este sistema de ecuaciones lineales:
26 64
"T R "N BR "F F R
37 75
26 64 +
=
d
1+ d
1 +
1
d
+1
+
+
1
b
1+
b
1 1 + +
37 26 75 64
1 + b
vd vs vB
37 75
Para identi…car completamente el modelo, se requiere, de restricciones adicionales al modelo. Estas restricciones estan dadas por el supuesto de como el banco central implementa su política monetaria. Bernanke y Mihov consideran los siguientes casos: 1. El banco central …ja la tasa de interés , esto es choques a la tasa de interés sólo re‡ejan choques de política monetaria: d = 1 y b = 1 2. El banco central …ja el nivel de las reservas no prestadas: d = 0 y b = 0
El modelo se puede estimar bajo cualquiera de estas dos especi…caciones para el esquema de política monetaria. En ambos casos, el regímen de política monetaria impone dos restricciones adicionales, por lo el modelo esta sobreidenti…cado, existe una restricción adicional en cada caso.Para estimar el modelo estructural Bernanke y Mihov siguen el siguiente procedimiento: 1. Estiman el sistema en forma reducida utilizando de MCO 2. Recuperan el vector de choques de la forma reducida 3. Estiman los parámetros de la forma estructural utilizando el método generalizado de momentos. El estimador de MGM se puede obtener minimizando la siguiende forma cuadrática:
bh
d = vec
i bh
B01 B01 vec
B01 B01
i
El modelo es su…cientemente ‡exible para permitir estimar y probar la existencia de cambios de regímenes monetarios, y recuperar una secuencia de choques de política monetaria consistente para todo el periodo muestral.
8.2 RESTRICCIONES CONTEMPORÁNEAS
95
Puesto que existen más momentos que parámetros a estimar, se pueden realizar una prueba de hipótesis sobre la validez de los momentos. Si los momentos son adecuados, d 0, si este no es el caso, el valor de d sería mayor a cero. Bajo la hipótesis nula de validez de instrumentos d 2 (q ) donde, q representa el número de momentos.
8.2.2. Cushman y Zha Estos dos autores estiman un VAR estructural para medir los efectos de política monetaria en una economía pequeña y abierta. Usan información mensual de Canada correspondiente al periodo 1974-1993. El vector de variables endógenas que utilizan en el model esta dado por: Y t = [y ; P ; R ; W P X ; T m ; T x ; Y ; P ; R ; M1; exch]
donde:
y : Indice de producción industrial desestacionalizado de los EE.UU P : Índice de precios al consumidor de los EE.UU. R : Tasa de interés de los fondos federales, EE.UU W P X : Indice de precios de exportación mundia. T m : importaciones Canadienses de los EE.UU T x : Exportaciones Canadienses a los EE.UU Y : Índice de producción industrial de Canada P : Índice de precios al consumidor de Canada R : La tasa de interés de los bonos del tesoro de Canada a 3 meses M 1 : Agregado monetario M1 de Canada exch : tipo de cambio dólar canadiense versus dólar americano.
El modelo en forma reducida esta dado por:
donde Y t =
24 35 Y 1;t Y 2;t
A(L)Y t = " t
, A(L) =
24
A11 (L) A12 (L) 0
A22 (L)
35
Esto es, Y 1;t corresponden a las variables externas, y Y 2;t a las variables domésticas. El supuesto de matriz triangular superior de A(L) implica que las variables externas son exógenas, esto es choques a la economía Canadiense no tienen efectos en la economía de los EE.UU. Notese que estas restricciones no se usan para identi…car el modelo. Para identi…car los parámetros de la forma estructural se imponen una estructura recursiva sobre la matriz de relaciones contemporáneas entre variables externas y domésticas. Estos es, se asume que las variables externas no reaccionan a choques en variables domésticas, y tambien que el producto,
96
CAPíTULO 8 VAR ESTRUCTURALES
la in‡ación, las exportaciones e importaciones estan relacionadas contemporaneamente de forma recursiva. Adicionalmente, se consideran las siguientes restricciones estructurales: 1. Una ecuación de demanda por dinero: 2. Una curva de oferta de dinero 3. Un sector de información
8.3. Variables no estacionarias Cuando el conjunto de variables bajo análisis no son estacionarias, estas deben previamente transformarse en estacionarias tomando primeras diferencias. Supongamos el vector Y t contiene series no estacionarias y que Y t , transforma el vector de series a estacionarias. Las correspondientes representaciones en forma reducida y estructural del model VAR estan dados por: (L)Y t = o + t (L)Y t = o + V t
Donde: o = B01 o (L) = B01 (L) = B01 B01
´Donde, es una matriz diagonal En este caso, la representación de medias móviles para la transformación estacionaria esta dado : Y t = o + (L) t Y t = o + (L)V t
Donde:
97
8.3 VARIABLES NO ESTACIONARIAS (L) = (L)B01
Sin embargo, dado que las series de tiempo son no estacionarias, también existe una representación de medias móviles para las series en niveles. Beverdge y Nelson (1981) demuestran que cualquier serie de tiempo no estacionaria puede escribirse como la suma de dos componentes no observables, un componente que representa la tendencia estocástica y el otro el componente transitorio Y t = Y tP + Y tT : Stock and Watson (1988) and King Ploser Stock and Watson (1991) generalizan este resultado al caso multivariante. En este caso, Y t puede representarse como:
Y t = + Art + D (L) "t
(8.1)
Donde el componente permanente esta dado por Y tP = + Art , donde rt esta de…nido como un paseo aletorio: rt = + rt1 + t
y el componente transitorio esta dado por Y tT = D (L) "t . KPSW (1991) muestran que la representación de tendencias comunes puede derivarse de la representación de medias móviles del vector Y t : Si Y t tiene un representación de Wold dada por Y t = + (L)"t ;
se puede resolver iterativamente para Y t Y t = Y 0 + o t + (1) t + t
0
(8.2)
P
donde t y t se de…nen como: t = ("1 + "2 + "3:: + "t ) y t = s "t+s donde s = (s+1 + s+2 + s+3 ) es una secuencia …nita. Note que (1) mide el impacto de largo plazo de los diferentes choques , cuando las series son estacionarias, (1) = 0, por tanto todos los choques solo tienen efectos transitorios sobre el nivel de la serie.
8.3.1. Restricciones de Largo Plazo Blanchard y Quah(1989) muestran que se pueden utilizar restricciones de largo plazo para identi…car modelos estructurales. Ellos utilizan un modelo bivariado para el producto y la tasa de desempleo que les permite identi…car
98
CAPíTULO 8 VAR ESTRUCTURALES
choques estructurales de demanda y de oferta. El modelo que utilizan se puede escribir como:
24 35 24 35 24 24 35 24 35 24 b11 b12 b21 b22 111
112
121 122
1
Y t
1o
=
1o
ut
Y t1 ut1
+
35
v1;t v2;t
+
35
Para identi…car el modelo requerimos estimar los cuatro parámetros bij Para ello podemos utilizar la información de la estimación de la versión en forma reducida del VAR. En particular sabemos que:
b b
= B01 B01
b
donde corresponde a la matriz de varianza covarianza de los residuos de la forma reducida y , la matriz de varianza covarianza de los residuos de la forma estructural. Suponiendo que = I
b
b
b
= B01 B01
Es un sistema de ecuaciones no lineales que permiten identi…car n(n+1) 2 parámetros, en el caso del modelo de BQ, esto nos dan 3 condiciones para estimar 4 parámetros. Por tanto necesitamos una restricción adicional. Esta restricción adicional se obtiene a partir de los efectos de largo plazo de los choques. La representación de medias móviles del modelo esta dada por: yt = 1 + 11 (L)v1;t + 12 (L)v2;t ut = 1 + 21 (L)v1;t + 22 (L)v2;t
Donde, v1;t representan choques de oferta y v2;t choques de demanda. Para identi…car el modelo BQ hacen los siguientes supuestos sobre el impacto de largo plazo de los choques: 1. Choques transitorios ( de demanda) no tienen efecto sobre el nivel de producto y desempleo. 2. Los choques de oferta si tienen efecto permanente en el producto pero no en el nivel de desempleo.
99
8.3 VARIABLES NO ESTACIONARIAS
La restricción de que choques de demanda no tienen efectos de largo plazo en el nivel de producto equivale a asumir que: 12 (1) =
1
X
s12 = 0
s=0
La restricción que los choques " 1;t y " 2;t no tienen efectos de largo plazo en u t equivale a: @u t+s = lm s2 j = 0 s!1 @v j;t s!1 lm
que se cumple automáticamente puesto que la tasa de desempleo es estacionaria. Las restricciones de BQ implican que la matrix de impacto de largo plazo es triangular inferior. (1) =
24
11 (1)
0
21 (1) 22 (1)
35
De la representación de medias móviles del modelo, tenemos que: Y t = (1) "t = (1) B01 V t = (1) V t
por lo tanto, debe cumplirse que : (1) B01 = (1)
Escribiendo esta última ecuación en forma extensiva tenemos:
24
11 (1) 12 (1) 21 (1) 22 (1)
35 24 35 24 b11 b12
=
b21 b22
11 (1)
0
21 (1) 22 (1)
35
De donde obtenemos la restricción adicional que permite identi…car todos los parámetros del modelo estructural. 11 (1)b11 + 12 (1)b12 = 0
Resumiendo, el modelo contiene cuatro condiciones para determinar cuatro parámetros del modelo estructural en función de los parámetros del modelo en forma reducida. b11 =
p
11
100
CAPíTULO 8 VAR ESTRUCTURALES
b21 =
b22 =
p 21
11
s p 22
21 11
11 (1)b11 + 12 (1)b12 = 0
Una vez identi…cado estos parámetros, la matriz B0 se utiliza para recuperar todos los parámetros estructurales, y la secuencia de choques de oferta y demanda. Conociendo la secuencia de choques de oferta y demanda se puede estimar el componente permanente del producto, y el componente transitorio, la tendencia y la brecha producto. Evaluando Neutralidad en el Largo plazo King y Watson (1997) resumen varios trabajos sobre VAR estructurales bivariados que buscan probar la hipótesis de neutralidad en el largo plazo. La neutralidad establece que cambios en variables nominales no tienen efectos reales en el largo plazo sobre variables reales. Algunas de las proposiciones de neutralidad en el largo plazo pueden son: 1. Un cambio permanente in la cantidad de dinero no tiene efectos de largo plazo sobre el nivel de producto real 2. Un cambio permanente en la tasa de in‡ación no tienen efectos en el largo plazo en el nivel de desempleo ( una curva de Phillips vertical). 3. Un cambio permanente en la tasa de in‡ación no tiene efectos permanentes en las tasas de interés real. ( Ecuación de Fisher) Para probar la neutralidad, se requiere de variables no estacionarias I (1) . Por ejemplo si queremos probar al primera de las hipótesis anteriores, podemos establecer: yt = 1 + 11 (L)"1;t + 12 (L)"2;t mt = 1 + 21 (L)"1;t + 22 (L)"2;t
Donde, "1;t representan choques exogenos al producto que no estan correlacionados con los choques exogenos a la cantidad de dinero, "2;t : Para
101
8.3 VARIABLES NO ESTACIONARIAS
probar la hipótesis de si el dinero es neutral, debemos responder la pregunta, ¿ llevan cambios exogenos en la cantidad de dinero "2;t a cambios permanentes en yt ? De la descomposición de Beverdge y Nelson sabemos que el efecto de largo plazo de las varibles esta dada por la siguiente matriz
24 35 24 1;t 2;t
=
11 (1) 12 (1) 21 (L) 22 (1)
35 24 X 35 X "1;t "2;t
donde, 1;t y 2;t representa las tendencias estocásticas que afectan al producto y la cantidad de dinero. La hipótesis de nuestralidad equivale a probar que: 12 (1) = 0, esto es los choques monetarios solo tienen efectos transitorios en el producto. La tendencia del producto esta únicamente determinada por choques reales, por ejemplo choques de productividad, " 1;t Para realizar esta prueba de hipótesis primero tenemos estimar los parámetros de la representación de medias móviles. Para ello debemos seguir los siguientes pasos: 1. Estimar los parámetros de la representación VAR. 2. Invertir a partir de los parámetros de la representación VAR, los parámetros de medias móviles
Capítulo 9
Modelos VAR en la práctica 9.1. Introducción En este apartado desarrollaremos en forma empírica la estimación de un sistema de Vectores autorregresivos (VAR), enfatizando el análisis intuitivo de los resultados y la rigurosidad de la construcción de los procedimientos. Para ello utilizaremos el entorno Eviews 5.1., que nos permitirá trabajar el tema con suma facilidad. Al …nal del capitulo se expone la estimación de estos modelos en MATLAB.
9.2. Motivación Muchas veces cuando queremos realizar un análisis riguroso de Politica Económica de tipo macro llegamos a la conclusión de que es necesario un modelo estructural que permita especi…car relaciones dinámicas entre varias variables. La implementación de dicho modelo nos debe proporcionar información sobre como reaccionaría una variable macroeconómica (exógena) ante un choque de politica macroeconómica (endógena). El interés deberá radicar más que nada en cuál es la dirección del efecto de este choque ó respuesta, de modo que se encuentre acorde con la teoría, es decir, que los resultados numéricos repliquen la intuición derivada de la teoría y no los contradigan. De otra forma, será necesario revisar la especi…cación del modelo. A lo largo del siglo XX el desarrollo de la macroeconomía presentó diversas variantes, donde una de las más importantes fue la Crítica de Lucas (1976), donde se señalaba principalmente la invalidez de los modelos macroeconómicos utilizados hasta ese entonces, dado que no re‡ejaban las rela103
104
CAPíTULO 9 MODELOS VAR EN LA PRÁCTICA
ciones macroeconómicas correctas, no tenían poder predictivo y por tanto no eran útiles para efectos de evaluación de política económica ( policy ). El principal problema radicaba en la especi…cación de estos modelos, donde no se observaban relaciones dinámicas sino mas bien estáticas. Es decir, toda relación entre variables era principalmente de tipo contemporáneo, con parámetros estimados constantes en el tiempo. En este escenario es fácil entender por qué toda evaluación de política era susceptible de ser errada, debido a que las relaciones especi…cadas distaban signi…cativamente de las correctas y por ende las predicciones y los escenarios futuros eran completamente inciertos. La crítica se basa en que la econometria clásica asumía parámetros estimados como reales y es más, constantes en el tiempo. Además, dichos modelos no incorporaban el componente de expectativas macroeconómicas (valores esperados) como determinante potencial de valores futuros en variables macroeconómicas, no expresaban relaciones dinámicas (intertemporales) y presentaban parámetros inestables. En suma, fue necesario desarrollar una nueva forma de hacer econometría, de forma tal que permita realizar evaluaciones de política, predecir escenarios y evaluar relaciones dinámicas. Una de estas alternativas fue el sistema de vectores autorregresivos discutido en primera instancia por Sims (1980) 1 y que desarrollamos en este capítulo como una alternativa plausible para identi…car el mecanismo de transmisión de política monetaria.
9.3. El sistema de Vectores autorregresivos: Una Recapitulación Un modelo VAR es un sistema de ecuaciones de n variables macroeconómicas relacionadas entre si. La idea es identi…car las respuestas de cada una de las variables frente a variaciones en las n 1 restantes. En modelos más avanzados se agrupan las variables entre las que son de política (que van a presentas perturbaciones) y aquellas que no son de política (y que no interesa saber su respuesta ante la perturbación anterior. Ello se verá más adelante.
9.3.1. Estimación Comunmente el sistema VAR se expresa de la siguiente forma: 1
Sims (1980): "Macroeconomics and Reality", Econometrica, Vol. 48, No. 1 (Jan., 1980) , pp. 1-48
9.3 EL SISTEMA DE VECTORES AUTORREGRESIVOS: UNA RECAPITULACIÓN 105
26 66 4
y1t y2t ::: ynt
37 77 5
=
26 64
A11 (l) ::: A1n (l) :::
:::
:::
An1 (l) ::: Ann (l)
37 26 75 664 X
y1t y2t ::: ynt
p
donde Aij (l) toma la forma de
37 77 5
+
26 66 4
C 1 C 2 ::: C n
37 77 5
+
26 66 4
"1t "2t ::: "nt
37 77 5
aijk lk , donde l es el operador de rezago
k=1
y p es el número de rezagos …jado por el modelador. Este sistema expresa un set de relaciones entre los valores presentes de cada variable con los valores pasados de si misma y de las restantes. Estos modelos son atractivo porque permiten expresar en forma sencilla relaciones intertemporales entre variales macroeconómicas y en base a ello analizar la secuencialidad de los mecanismos de transmisión de política, así como la causalidad entre variables (causalidad à la Granger ). En forma matricial el modelo se expresa como:
Y t = o + 1 Y t1 + 2 Y t2 :::: + p Y t p + t
con t N (0; ) Luego el modelo anterior se expresa en forma compacta apilando las variables exógenas, en este caso predeterminadas en un vector de variables X t tal que:
X t =
26 66 4
1 Y t1
.. .
Y t p
37 77 5
y = [o 1 ::: p ]
entonces, el VAR(p) puede escribirse como: Y t = 0 X t + t
En el capitulo anterior se demostró que la estimación Máximo Verosímil de este modelo coincidía con la regresión MCO de cada ecuación, es decir, esta última es consistente y de mínima varianza. Por ello, tomaremos este hecho como dado y procederemos a la estimación bajo MCO. El estimador MCO viene dado por:
106
CAPíTULO 9 MODELOS VAR EN LA PRÁCTICA
"X #"X # "X #"X # T
0
k (kp+1) =
T
Y t X 0 t
t=1
Asimismo se tiene que:
X t X 0 t
t=1 t
T
j10(kp+1) =
T
y j;t X t0
t=1
1
X t X t0
1
t=1 t
9.3.2. Función de impulso-respuesta Sea un modelo V AR( p) en forma reducida como el de…nido líneas arriba, la descomposición de Wold vendrá dada por: Y t = 0 + t + 1 t1 + 2 t2 + :::::
Sea entonces B0 la matriz que pre-multiplica a Y t en la forma estructural y V t los errores de la misma, la representac‘ón en medias móviles se convierte en:
Y t = c + B01 V t + 1 B01 V t1 + 2 B01 V t1 + :::
Por tanto, la función de impulso-respuesta describe la reacción del vector de variables Y t+s ante un choque de una sola vez en la variable y jt , ceteris paribus . En estricto se mide como: @Y t+s = s b j @v j:t
Donde:
h i
B01 = b1 b2 : : : bk
9.3.3. Descomposición de Varianza
En estricto, este concepto señala la participación o contribución al Error de Predicción de cada uno de los n choques presentes en el VAR. La varianza del error de predicción puede escribirse como la suma de la contribución marginal a la varianza total de la varianza de cada choque:
9.4 EL MODELO DE STOCK Y WATSON (2001) 6 j=k
=
X
j=1
var(v j;t ) b j b j 0 + 1 b j b j 0 01 + :::s1 b j b j 0 0s1
107
Entonces, la contribución marginal del choque j en la varianza del error de predicción de la variable i esta dado por: 100
var(v j;t)diag
b j b j 0 + 1 b j b j 0 01 + :::s1 b j b j 0 0s1 (i; i)
A medida que el horizonte de predicción se hace mayor, la varianza del error de predicción converge a la varianza no condicional de la serie 0 de manera análoga a lo que demostramos en la clase de predicción para el caso de modelos univariantes.
9.4. El Modelo de Stock y Watson (2001)2 Este modelo consta de tres variables: in‡ación, tasa de desempleo y la tasa de interés de la Reserva Federal (FED). La primera especi…cación que haremos será la de un VAR no restringido (à la Cholesky), ordenando las variables desde la más endógena hasta la más exógena, es decir, la matriz que pre-multiplica al vector de variables endógenas es triangular superior. En este caso, basándonos en un esquema tradicional de política monetaria, el mecanismo de transmisión va desde la tasa de interés, pasando por la demanda agregada 3 y terminando en la variación del índice de precios 45 . Sean entonces: t : In‡ación (Variación porcentual del Índice de Precios) t : Tasa de Desempleo Rt : Tasa de Interés de la Reserva Federal (FED) 2
Stock y Watson (2001), “Vector Autorregressions”, Journal of Economic Perspectives (JEP), Volumen 15, Número 4, Páginas 101-115. 3 Existen muchas formas de representar la actividad económica, que pueden ir desde el mismo producto, una brecha estimada o bien la tasa de desempleo, que es la que se utiliza en este caso. Se deja abierta la posibilidad a incluir otra variable plausible. 4 Para mayores detalles ver: Bernanke & Blinder (1992), “The Federal Funds Rate and the Channels of Monetary Transmission”. American Economic Review, Vol.82, N o 4, 901921. 5 Luego se realizará el análisis de causalidad à la Granger.
108
CAPíTULO 9 MODELOS VAR EN LA PRÁCTICA
Tendremos que: Y t =
26 64
t t Rt
37 75
La especi…cación reducida del modelo V ar( p) es como sigue: Y t = A0 + A1 Y t1 + A2 Y t2 + ::: + A p Y t p + t
Donde cada matriz A i es de orden (nxn), siendo en este caso particular n = 3 y donde cada componente de t es una combinación lineal de los errores estructurales. Los datos disponibles son trimestrales y van desde 1947 (1T) hasta 2001 (3T), pero para obtener una muestra balanceada se tomarán solamente desde 1960 (1T) hasta 2000 (4T).
9.5. Var en Eviews Eviews cuenta con rutinas bastante amigables y especializadas en la estimación de los VARs. Por ello en la medida de lo posible trataremos de sacar el máximo provecho a dicha ventaja, centrándonos en la correcta especi…cación de los modelos, así como en el análisis de los resultados.
9.5.1. Estimación Contando con el archivo de datos (.wf1) es posible crear un programa que lo llame directamente y opere la estimación del VAR. A continuación ejecutaremos esta tarea:
’Especi…camos la ruta del work…le cd c:n ’Llamamos al work…le load stock_watson_v2
Con el work…le abierto es posible estimar el VAR reducido, no olvidando que es necesario encontrar el número óptimo de rezagos antes de llevar a cabo esta tarea. El criterio que utilizaremos será el de Schwarz, que entre varios más penaliza la inclusión de un rezago adicional en el VAR, siendo posible identi…car un modelo parsimonioso. En la práctica estimamos varios modelos VAR en forma consecutiva, de modo que aquel que tenga un Schwarz mayor al aumentar un rezago deberá ser eliminado, quedándonos con el inmediato anterior. ’Estimamos un VAR no restringido
9.5 VAR EN EVIEWS
109
!lagmax=6 vector(!lagmax) schw var var1.ls 1 1 in‡ation un_rate fed_funds @ c schw(1)=var1.@sc !i=2 var var1.ls 1 !i in‡ation un_rate fed_funds @ c schw(!i)=var1.@sc while schw(!i)>schw(!i-1) var var1.ls 1 !i in‡ation un_rate fed_funds @ c schw(!i)=var1.@sc !i=!i+1 wend d schw scalar lagopt=!i
Luego de estos comandos Eviews ha encontrado el número óptimo de rezagos y posteriormente lo ha almacenado en el objeto escalar (scalar) llamado lagopt. A continuación, dado el número de rezagos efectuamos el test de causalidad à la Granger. El resultado se muestra al …nal. ’Análisis de Causalidad à la Granger group g1 in‡ation un_rate fed_funds freeze(granger_test) g1.cause(!i) show granger_test
110
CAPíTULO 9 MODELOS VAR EN LA PRÁCTICA
A priori podemos asumir que las tres variables son efectivas para predecirse entre sí mismas, por lo que recurrimos al esquema teórico del mecanismo de transmisión monetario para tomar una decisión …nal. Adicionalmente hacemos que los resultados de esta selección de rezagos se muestren en una tabla, de tal forma que no queden dudas al respecto. Es más, hacemos que la tabla salga con un rezago adicional. !j=!i+1 freeze(lag_choice) var1.laglen(!j) show lag_choice
Los resultados se muestran en la siguiente tabla:
Los resultados del VAR no restringido se presentan a continuación: freeze(var_chol) var1.results show var_chol
111
9.5 VAR EN EVIEWS
Como sabemos, las estimaciones son realizadas por el método de Mínimos Cuadrados Ordinarios (MCO) y por ende los análisis de signi…cancia y residuos son análogos al análisis univariante, por lo que no haremos mucho hincapié al respecto. Cuando se habla de V AR el producto …nal más importante está formado por dos componentes: La función de impulso-respuesta y la Descomposición de Varianza. En el caso del primero, es necesario obtener la descomposición de Wold del modelo estimado, es decir, la representación V MA (1).
9.5.2. Función de impulso-respuesta En Eviews generamos la función de impulso-respuesta y gra…camos la reacción hasta por 24 períodos. Le añadimos bandas sumando y restando el doble de la desviación estándar. Recordando que:
@Y t+s @v j:t
= s b j
freeze(chol_imp) var1.impulse(24, se=a) show chol_imp
112
CAPíTULO 9 MODELOS VAR EN LA PRÁCTICA
Observamos detenidamente los resultados, recordando el ordenamiento de las variables ( t , t , Rt ). La primera columna muestra la respuesta de las tres variables (en ese orden) ante un choque en la in‡ación, la segunda ante un choque en la tasa de desempleo y la tercera ante un choque en las tasas de interés. Los signos deben ser los correctos, de tal forma que podamos asegurar que los datos describen correctamente el comportamiento de la economía agregada. Claramente observamos los trade-o¤ entre in‡ación y desempleo y entre tasa de interés y desempleo ( proxy de la actividad económica) en el corto plazo. De igual forma apreciamos la relación directa entre in‡ación y tasa de interés, etc. Para el tema de la signi…cancia de estos choques debe tenerse en cuenta lo siguiente: En general un choque signi…cativo presenta las dos bandas de con…anza en el mismo cuadrante, es decir, las dos del lado positivo o viceversa.
9.5.3. Descomposición de Varianza La contribución marginal de cada choque a la varianza total del vector de series endogenas esta dado por:
100 E (Y t+s
Y t+s=t)(Y t+s Y t+s=t) 1 var(v j;t)
b j b j 0 + 1 b j b j 0 01 + :::s1 b j b j 0 0s1
En Eviews podemos construir la Descomposición de Varianza tal como sigue:
f reeze(chol_decvar) var1.decomp(24,se=a) in‡ation un_rate fed_funds
9.5 VAR EN EVIEWS
113
show chol_decvar
El resultado nos muestra evolución de la contribución de cada variable a lo largo del horizonte de predicción.
La primera …la nos muestra la contribución de cada una de las tres variables (según el orden señalado antes) al error de predicción de la in‡ación. La segunda es para la predicción de la tasa de desempleo y la tercera para la tasa de interés. Observamos claramente que en los primero períodos solo es relevante el choque de la misma variable, lo que va yendo en detrimento a medida que el horizonte se hace más amplio. Es decir, en un primer momento funciona el modelo adaptativo y a medida que avanza el tiempo la construcción del modelo (estructura) comienza a tomar mayor ponderación.
9.5.4. Extensión: Predicción en VAR El proceso de predicción en un modelo VAR es una extensión al univariante, teniendo en cuenta la regla de formación para el horizonte dado. Sea el modelo V AR siguiente: Y t = A 0 + A1 Y t1 + A2 Y t2 + ::: + A p Y t p + t
La predicción para el siguiente período será: E t Y t+1 = A0 + A1 Y t + A2 Y t1 + ::: + A p Y t p+1
114
CAPíTULO 9 MODELOS VAR EN LA PRÁCTICA
Lo primero será parametrizar el programa, haciendo variar la muestra cuando cambia el horizonte a predecir. En esta oportunidad haremos la predicción de un período hacia delante (h = 1). En clases posteriores extenderemos este tema.
’Horizonte de Predicción !s=1 smpl 1960:01 2000:4-!s
Luego, en Eviews es necesario extraer las matrices A i para llevar a cabo esta operación: matrix pi=var1.@coefmat for !j=1 to 3 vector(3) _r{!j} next for !k=1 to !i matrix(3,3) pi{!k}t _r1=@rowextract(pi,!k) rowplace(pi{!k}t,_r1,1) _r2=@rowextract(pi,!k+!i) rowplace(pi{!k}t,_r2,2) _r3=@rowextract(pi,!k+2*!i) rowplace(pi{!k}t,_r3,3) matrix (3,3) pi{!k} pi{!k}=@transpose(pi{!k}t) d pi{!k}t next vector(3) mut mut=@rowextract(pi,!k+2*!i) vector(3) mu mu=@transpose(mut) d _r* mut
Luego convertimos las variables a matriz y operamos la predicción: scalar _N=var1.@regobs !n=_N matrix _Yt=@convert(g1) matrix _Y=@transpose(_Yt) d _Yt for !j=0 to !i-1 !f=!j+1 vector _Y{!f}=@columnextract(_Y,!n-!j) next
9.6 VAR EN MATLAB 8
115
vector(3) _Yf{!s} _Yf{!s}=mu+pi1*_Y1+pi2*_Y2
9.6. VAR en MATLAB7 9.6.1. Estimación Recordando el modelo de Stock y Watson procedemos a replicarlo en el entorno MATLAB. La idea está en estimar por mínimos cuadrados ordinarios (MCO) los parámetros del modelo siguiente, para luego hacer el análisis de impulso-respuesta y descomposición de varianza. Y t = A 0 + A1 Y t1 + A2 Y t2 + ::: + A p Y t p + t
Ya tenemos nuestro primer objetivo, para lo cual debemos implementar una rutina que replique dicho procedimiento en forma ordenada y plausible. Antes de empezar a programar debemos tener bien claro el objetivo ó lo que queremos lograr, de tal forma que no perdamos tiempo de trabajo. En primer lugar necesitamos los datos para poder estimar el sistema V AR, los cuales reán importados directamente de una hoja Excel. Para ello se debe …jar en dicho arhivo las variables a utilizar, teniendo cuidado en seguir el orden utilizado en la estimación por Eviews. Una vez que se tenga este archivo listo se debe ir a la opción ’FilenImport data...’, para luego dar click en ’…nish’. Con esta sentencia MATLAB creará una matriz llamada ’data’ que contendrá las series de tiempo para poder ser trabajadas. Luego de esto, esposible pasar al paso de estimación, capturando los datos importados: y=data;
En resumidas cuentas necesitamos una función que muestre los resultados de la estimación del sistema V AR, contando como argumentos principales a los datos yt , el número de rezagos p y si es necesario un grupo de variables exógenas, tema que se verá más adelante. Por simplicidad asumiremos que para este caso no existen variables exógenas, dejando abierta la posibilidad para temas futuros de análisis. 7
La base de estos programas puede ser encontrada en "Econometric Toolbox for MATLAB"en http://www.spatial-econometrics.com
116
CAPíTULO 9 MODELOS VAR EN LA PRÁCTICA
Entonces, concluimos que necesitamos una función de dos argumentos tal como: function results = vare(y,p)
La primera línea del programa está escrita, con lo cual ya tenemos el ob jetivo trazado. En este caso y representa la matriz principal de datos (todas las series de tiempo en columnas concatenadas en una sola matriz). Luego, el resultado ó imagen de la función deberá tomar una forma particular, dado que tenemos varias cosas que exponer cuando estimamos un var, por ejemplo, el criterio de información Schwarz, así como la signi…cancia global, etc. En suma, es necesario construir un ob jeto que cumpla con las propiedades de almacenar varios resultados a la vez, así como la descripción de los mismos. En primera instancia el modelador podria escoger una matriz cono el objeto deseado, pero dado que sería mejor incorporar la descripción de los resultados, es conveniente trabajar con texto y números a la vez. Para ello existe el objeto estructuras ( struct), que es capaz de crear una matriz ó arreglo de dos columnas distintas: una para la descripcion (texto) y la segunda para los resultados númericos ó strings (texto). Es decir, el elemento struct es una suerte de tabla de resultados …jada a gusto del modelador. En realidad el término estructura es mucho más complejo que una tabla o una matriz, pero por simplicidad estos términos nos facilitan el entendimiento del procedimiento. Luego de obtener los resultados, estos podrán ser llamados como se hace comunmente con las componentes de un vector, tal como results(i), suponiendo que en la posición i se encuentra guardado el resultado tal como, por ejemplo, criterio de Schwarz. No´tese que en MATLAB ya existe una función var para calcular la varianza de un proceso estocástico (vector), por lo que será necesario escoger un nombre distinto, tal como vare. La estimación de un sistema VAR está basada en la lógica de mínimos cuadrados ordinarios (MCO), por lo que será necesario crear una función que replique este procedimiento. Asimismo, existen varias funciones adicionales que seran utilizadas en el desarrollo de este procedimiento y que se encuentran descritas en el apéndice de este capítulo. Continuando con la función creada, capturamos el número de observaciones (T) y el número de ecuaciones a estimar (número de variables del sistema var). [T n] = size(y);
Dentro de la estructura se especi…ca el nombre de la función utilizada en la estimación, que siempre es útil cuando se trabaja con diversos experimentos computacionales simultáneamente. results.meth = ’vare’;
9.6 VAR EN MATLAB 10
117
Dado el número de rezagos propuesto, obtenemos el número de observaciones disponibles para la estimación ( Te). Asimismo, se agrega el término constante como variable explicativa: Te = T - p; k = n*p + 1; nvar = k; results.nvar = nvar;
Como siempre es necesario tener a los rezagos de las variables en forma matricial, la función mlag permite obtener este tpo de matriz, que facilita enormemente los cálculos 9 . xlag = mlag(y,p);
Guardamos la información necesaria en la estructura de resultados: results.T = Te; results.n = n; results.p = p;
Construimos la matriz X t , que representa las variables explicativas del modelo en cuestión: xmat = [xlag(p+1:T,:) ones(T-p,1)];
Hasta aqui hemos recolectado la información relevante necesaria para la estimación de un VAR por MCO. El siguiente paso será crear un secuencia iterativa que estime los parámetros de cada ecuación por MCO y que construya los indicadores relevantes para el diagnóstico de la regresión, tales como R2 , el estadístico F , el criterio de Schwarz, etc. El objeto results es una estructura que permite realizar dicho proceso sin mayores di…cultades: Dado que es un proceso iterativo, especi…caremos el código para estimar una sola ecuación y luego este será repetido n veces. for j=1:n;
En el ejercicio, la variable j es la endógena en la estimación por MCO: yvec = y(p+1:T,j);
Estimando por MCO la ecuación correspondiente: res = ols(yvec,xmat);
La función ols construye una estructura de resultados ( res) similar a la que deseamos obtener en este ejercicio, cuyo detalle se encuentra en el apéndice a este capítulo. En este caso capturamos los parámetros estimados y los valores de los estadísticos t correspondientes para la ecuación j y los almacenamos en nuestra estructura de resultados: results(j).beta = res.beta; results(j).tstat = res.tstat; 9
Ver el detalle de esta y otras funciones en el apéndice de este capítulo.
118
CAPíTULO 9 MODELOS VAR EN LA PRÁCTICA PRÁCTICA
De igual forma, construimos las probabilidades asociadas a estos valores y los registramos en la estructura de resultados: tstat = zeros(nvar,1); tstat = res.tstat; tout = tdis_prb(tstat,Te-nvar); results(j).tprob = tout;
Transferimos el resto de indicadores de la regresión MCO realizada: Residuos: results(j).resid = res.resid;
Suma de cuadrados de los residuos: sigu = res.resid’*res.resid;
Valores estimados: results(j).yhat = res.yhat;
Valores observados: results(j).y = yvec;
R-cuadrado: results(j).rsqr = res.rsqr; res.rsqr; % r-cuadrado
R-cuadrado ajustado: results(j).rbar = res.rbar;
Varianza estimada: results(j).sige = res.sige;
El estadístico Box-Pierce Q para los residuos: elag = mlag(res.resid,p); etrunc = elag(p+1:Te,:); elag(p+1:Te,:); rtrunc = res.resid(p+1:Te,1); res.resid(p+1:Te,1); qres = ols(rtrunc,etrunc);
Ajustando por el numero de rezagos: if p ~= 1 boxq = (qres.rsqr/(p-1))/((1-qres.rsqr)/(T (qres.rsqr/(p-1))/((1-qres.rsqr)/(Te-p)); e-p)); else boxq = (qres.rsqr/(p))/((1-qres.rsqr)/(T (qres.rsqr/(p))/((1-qres.rsqr)/(Te-p)); e-p)); end; results(j).boxq = boxq;
Construyendo los estadísticos de signi…cancia global: En este apartado se mide la signi…cancia por cada rezago, excluyendo uno por cada iteración: for r=1:n; xtmp = []; for s=1:n; if s ~= r xlag = mlag(y(:,s),p);
9.6 VAR EN MATLAB 11 xtmp = [xtmp trimr(xlag,p,0)]; end; end; xtmp = [xtmp ones(Te,1)]; b = xtmpnyvec;
etmp = yvec-xtmp*b; sigr = etmp’*etmp; ftest(r,1) = ((sigr - sigu)/p)/(sigu/(Te-k)); end; results(j).ftest = ftest; results(j).fprob results(j).fprob = fdis_prb(ftest,p,T fdis_prb(ftest,p,Te-k); e-k); end;% end; % Fin del programa programa
Mostramos los resultados: prt_var(results,vnames)
***** Vector Autoregressive Model ***** Dependent Variable = inflation R-squared = 0.8341 Rbar-squared Rbar-squared = 0.8277 sige = 1.0906 Q-statistic = 9.7705 Nobs, Nvars = 162, 7 *************************************** ******************* **************************************** *************************** ******* Variable Coefficient t-statistic t-probability inflation lag1 0.700602 8.862242 0.000000 inflati inflation on lag2 0.224005 0.224005 2.700673 2.700673 0.007690 0.007690 n_ra te te la g1 g1 - 0. 0.661919 - 1. 1. 95 950173 0.052960 n_rate lag2 0.568061 1.756862 0.080916 fed_fund lag1 lag1 0.130823 1.321686 0.188220 fed _f _fu nd nd la g2 g2 - 0. 0. 09 094 62 622 - 0. 0. 93 93 50 50 41 41 0. 35 35 12 122 2 constant 0.609946 1.710357 0.089201 ****** Granger Causality Tests ******* Variable F-value Probability inflation 189.084954 0.000000 n_rate 2.309166 0.102752 fed_fund 1.199544 0.304113 0.304113
119
120
CAPíTULO 9 MODELOS VAR EN LA PRÁCTICA PRÁCTICA Dependent Variab le = un_rate R-squared R-squared = 0.9780 Rbar-squared Rbar-squared = 0.9771 sige = 0.0526 Q-statistic = 0.0550 Nobs, Nvars = 162, 7 ************************************** ********************************************************** **************************** ******** Variable Coefficient t-statistic t-probability inflation inflation lag1 0.035764 2.060187 0.041053 infla ti tio n la g2 g2 - 0. 0.028532 - 1. 1.566536 0.119262 n_rate n_rate lag1 1.520327 1.520327 20.398499 20.398499 0.000000 0.000000 n_ra te te la g2 g2 - 0. 0. 57 575917 - 8. 8.111388 0.000000 fed_fund fed_fund lag1 lag1 0.011441 0.526398 0.599364 fed_fund fed_fund lag2 lag2 0.018993 0.854726 0.394022 constan constantt 0.0992 0.099242 42 1.2673 1.267314 14 0.2069 0.206944 44 ****** Granger Causality Tests ******* Variable F-value Probability inflation 2.143652 0.120689 n_rate 2435.118546 0.000000 fed_fund fed_fund 6.905810 0.001340 0.001340
Dependent Variab le = fed_fu nd R-squared R-squared = 0.9177 Rbar-squared Rbar-squared = 0.9145 sige = 0.8690 Q-statistic = 16.1060 Nobs, Nvars = 162, 7 ************************************** ********************************************************** **************************** ******** Variable Variable Coefficient Coefficient t-statisti t-statisticc t-probabil t-probability ity infla ti tio n la g1 g1 - 0. 0.019923 - 0. 0.282335 0.778063 inflation inflation lag2 0.193685 2.616006 0.009776 n_ra te te la g1 g1 - 1. 1. 35 355610 - 4. 4.474370 0. 00 000015 n_rate lag2 1.272148 4.407671 0.000019 fed_fund fed_fund lag1 lag1 0.920866 10.422457 0.000000 fe d_ d_ fu fun d la g2 g2 - 0. 0. 03 03 18 18 34 34 - 0. 0. 35 35 24 242 2 0 .7 .7 25 250 01 01 constant 0.567454 1.782607 0.076608 ****** Granger Causality Tests ******* Variable F-value Probability inflation 8.644371 0.000276 n_rate 10.010499 0.000082 0.000082 fed_fund ed_fund 381.41 381.41282 2820 0 0.0000 0.000000 00
9.6.2. Función de impulso-res impulso-respuesta puesta Calculamos las funciones de impulso respuesta de la descomposición de Cholesky: [out1,out2]=irf(results,10,’o1’,vnames)
9.6 VAR EN MATLAB 12
121
122
CAPíTULO 9 MODELOS VAR EN LA PRÁCTICA
Como hemos podido apreciar, los resultados son similares a los reportados por Eviews 5.1., lo cual valida la rutina utilizada anteriormente. Debe recordarse que lo importante es …jar el ordenamiento de las variables en el archivo de datos adjunto, dado que esto afecta la estimación e identi…cacion del V AR.
Capítulo 10
Modelos Var Estructurales en la Práctica 10.1. Introducción En este capítulo replicaremos dos de los ejemplos más populares en materia de Vectores autorregresivos estructurales. Ambos fueron presentados en el capítulo anterior, por lo que en este apartado entraremos directamente a la aplicacón. El primer modelo lo presentan Stock y Watson (2001), es decir, haremos la especi…cación estructural del modelo del capítulo anterior. El segundo modelo es un Var simple bi-variado que utiliza restricciones de largo plazo para la identi…cacion de modelos estructurales, que se desarrolla con …nes didácticos. Finalmente veremos a Bernanke y Mihov (1998), quienes plantean un modelo del mercado interbancario, con el objetivo de identi…car choques de política monetaria.La mayor parte de este apartado se encuentra basada en el capitulo 6 de Favero (2001). 1 En el desarrollo del capítulo se ofrecen rutinas de Eviews 5.1. para la estimacion de estos modelos. Asimismo, al …nal de cada modelo se presenta la estimación del mismo ejercicio por medio de rutinas de MATLAB.
10.2. El Modelo de Stock y Watson (2001) De acuerdo al modelo resuelto en el capítulo anterior, especi…camos una forma diferente a la recursiva, restringiendo las relaciones contemporáneas entre variables. En este caso, asumiremos que la in‡ación no reacciona ante 1
Favero, Carlo (2001): .A pplied Macroeconometrics". Oxford University Press.
123
124CAPíTULO 10 MODELOS VAR ESTRUCTURALES EN LA PRÁCTICA un choque en la tasa de interés en el mismo período, lo que llevará a que la forma estructural solamente cuente con 5 parámetros a estimar. El sistema quedará sobre-identi…cado, por lo que será necesario validar esta especi…cación mediante un test de sobre-identi…cación, conocido como el estadístico J2 . De esta forma tendremos:
A =
0B B@
1
0
0
NA
1
0
0
NA 1
1C CA
, B =
0B B@
NA
0
0
0
NA
0
0
0
NA
1C CA
De esta forma, el sistema tiene 5 parámetros estructurales frente a 6 de la forma reducida, por lo que cuenta con un grado de libertad. Fijamos las matrices en Eviews: ’Especi…cación estructural matrix(3,3) pata pata.…ll(by=r) 1,0,0, na,1,0, 0,na,1 matrix(3,3) patb patb.…ll(by=r) na,0,0, 0,na,0, 0,0,na pata(3,1)=0 ’La in‡ación no reacciona en forma contemporánea a la tasa de interés
El número óptimo de rezagos se encuentra grabado en el objeto var, dado que fue utilizado en la última estimación por MCO (forma reducida). Asimismo, utilizamos los residuos de la forma reducida como condiciones iniciales para el algoritmo: freeze(svar1) var1.svar(rtype=patsr,conv=1e-5,namea=pata,nameb=patb,F0=residuos) ’show svar1 freeze(svar1_imp) var1.impulse(24, se=a,imp=struct) show svar1_imp freeze(svar1_decvar) var1.decomp(24,se=mc,imp=struct) in‡ation un_rate fed_funds ’show svar1_decvar var1.makeresids sres_in‡ation sres_un_rate sres_fed_funds group sresiduos sres_in‡ation sres_un_rate sres_fed_funds freeze(struct_resid) var1.resids ’show struct_resid 2
Estadístico comúnmente utilizado en la estimación por el Método Generalizado de Momentos (MGM)
10.2 EL MODELO DE STOCK Y WATSON (2001)
125
El test Chi-cuadrado acepta la hipótesis nula de sobreidenti…cación, lo cual implica que la especi…cación es correcta y no será necesario descartarla. Es decir, podría entrar como un modelo candidato para ser escogido entre
126CAPíTULO 10 MODELOS VAR ESTRUCTURALES EN LA PRÁCTICA otros más.
Response to Structural One S.D. Innovations ± 2 S.E. Response of INFLATION to Shock1
Response of INFLATION to Shock2
Response of INFLATION to Shock3
1.2
1.2
1.2
0.8
0.8
0.8
0.4
0.4
0.4
0.0
0.0
0.0
-0.4
-0.4
-0.4
-0.8
-0.8 5
10
15
20
-0.8 5
Response of UN_RATE to Shock1
10
15
20
5
Response of UN_RATE to Shock2 .5
.5
.4
.4
.4
.3
.3
.3
.2
.2
.2
.1
.1
.1
.0
.0
.0
-.1
-.1
-.1
-.2
-.2
-.2
-.3
-.3 10
15
20
10
15
20
5
Response of FED_FUNDS to Shock2 1.2
1.2
0.8
0.8
0.8
0.4
0.4
0.4
0.0
0.0
0.0
-0.4
-0.4
-0.4
-0.8
-0.8
-0.8
-1.2
-1.2 10
15
20
10
15
20
Response of FED_FUNDS to Shock3
1.2
5
20
-.3 5
Response of FED_FUNDS to Shock1
15
Response of UN_RATE to Shock3
.5
5
10
-1.2 5
10
15
20
5
10
15
20
Los resultados en este caso señalan que la actividad económica si reacciona ante cambios en la tasa de interés, no así la in‡ación, que se explica por los rezagos de política. Por otro lado, la tasa de interés reacciona ante cambios en la in‡ación, lo que podría estar relacionado con una regla de Taylor que mira hacia atrás para …jar sus decisiones de política. Del mismo modo, esta tasa reacciona ante cambios en la actividad económica, que se puede explicar por la misma a…rmación anterior.
10.3 EL MODELO CON RESTRICCIONES DE LARGO PLAZO 127 Variance Decomposition Percent INFLATIONvariance due to Shock1
Percent INFLATIONvariance due to Shock2
PercentINFLATIONvariance due toShock3
100
100
100
80
80
80
60
60
60
40
40
40
20
20
20
0
0 2
4
6
8
10
12
14
16
18
20
22
24
0 2
Percent UN_RATE variance du eto Shock1
4
6
8
10
12
14
16
18
20
22
24
2
Percent UN_RATE variance due to Shock2 100
100
80
80
80
60
60
60
40
40
40
20
20
20
0 2
4
6
8
10
12
14
16
18
20
22
24
4
6
8
10
12
14
16
18
20
22
24
2
Percent FED_FUNDS variance due to Shock2 90
90
80
80
80
70
70
70
60
60
60
50
50
50
40
40
40
30
30
30
20
20
20
10
10
10
0
0 4
6
8
10
12
14
16
18
20
22
24
10
12
14
16
18
20
22
24
4
6
8
10
12
14
16
18
20
22
24
Percent FED_FUNDS variance due to Shock3
90
2
8
0 2
Percent FED_FUNDS variance due to Shock1
6
Percent UN_RATEvariance due toShock3
100
0
4
0 2
4
6
8
10
12
14
16
18
20
22
24
2
4
6
8
10
12
14
16
18
20
22
24
En el corto plazo las predicciones son en base a modelos adaptativos, pero en la medida en que el horizonte de predicción es mayor se apela al uso de modelos estructurales, empleando otras variables para llegar a un mejor resultado.
10.3. El modelo con restricciones de largo plazo Este modelo utiliza las restricciones de largo plazo como una alternativa para obtener un sistema VAR exactamente identi…cado. En esta oportunidad utilizaremos un modelo sencillo, donde solamente estará presente la brecha de producto y la tasa de desempleo para Estados Unidos, en el período 1948-1987 (datos trimestrales). La idea es ilustrar con este ejemplo sencillo la utilidad de las restricciones de largo plazo para identi…car un Var estructural. Si bien es cierto un modelo bi-variado no necesita otra especi…cación sino la triangular, no haremos uso
128CAPíTULO 10 MODELOS VAR ESTRUCTURALES EN LA PRÁCTICA de dicho supuesto, asumiendo en cambio que una no afecta en el largo plazo a la otra. Sea un modelo V AR como el siguiente:
Y t = A0 + A1 Y t1 + A2 Y t2 + ::: + A p Y t p + t
donde: Y t =
0@ 1A yt ut
, yt es el producto real en diferencias y ut es la
tasa de desempleo. Lo que queremos demostrar es que ‡uctuaciones de corto plazo en el desempleo (pequeñas desviaciones) no afectan al crecimiento del producto en el largo plazo. Estimamos un Var reducido similar al modelo de Stock y Watson: ’ Programa que aplica la metodología de Restricciones de Largo Plazo en un Var Estructural ’ Fernando Pérez Forero ’ 11/01/06 ’===================================================== ’Especi…camos la ruta del work…le cd e:n ’Llamamos al work…le load blanquah ’Estimamos un Var no restringido !lagmax=6 vector(!lagmax) schw smpl 1948:1 1987:4 var var1.ls 1 1 dy u @ c schw(1)=var1.@sc scalar _N=var1.@regobs !n=_N !i=2 var var1.ls 1 !i dy u @ c schw(!i)=var1.@sc while schw(!i) >schw(!i-1) var var1.ls 1 !i dy u @ c schw(!i)=var1.@sc !i=!i+1 wend
10.3 EL MODELO CON RESTRICCIONES DE LARGO PLAZO 129 d schw scalar lagopt=!i ’Análisis de Causalidad à la Granger group g1 dy u freeze(granger_test) g1.cause(!i) ’show granger_test !j=!i+1 freeze(lag_choice) var1.laglen(!j) ’show lag_choice freeze(chol_imp) var1.impulse(24, se=a) show chol_imp freeze(chol_decvar) var1.decomp(24,se=a) dy u ’show chol_decvar var1.makeresids res_dy res_u group residuos res_dy res_u freeze(chol_resid) var1.resids ’show chol_resid freeze(var_chol) var1.results ’show var_chol freeze(var_chol_eq) var1.representations
Recogemos los residuos de la forma reducida para utilizarlos como condiciones iniciales. Luego, especi…camos la forma estructural con restricciones de largo plazo y estimamos el modelo: ’Especi…cación estructural var1.cleartext(svar) var1.append(svar) @lr2(@u1)=0 freeze(results) var1.svar(rtype=text,conv=1e-4,F0=residuos) show results freeze(svar_imp) var1.impulse(24,se=a,imp=struct,a) dy u show svar_imp freeze(svar1_decvar) var1.decomp(24,se=a,imp=struct) dy u show svar1_decvar var1.makeresids sres_dy sres_u group sresiduos sres_dy sres_u freeze(struct_resid) var1.resids ’show struct_resid
Los resultados del modelo se presentan a continuación:
130CAPíTULO 10 MODELOS VAR ESTRUCTURALES EN LA PRÁCTICA
10.4 MIDIENDO LOS CHOQUES DE POLíTICA MONETARIA: BERNANKE Y MIHOV (1 Accumulated Response to Structural One S.D. Innovations Accumulated Response of DY to Shock1
Accumulated Response of DY to Shock2
.008
.008
.004
.004
.000
.000
-.004
-.004
-.008
-.008
-.012
-.012 2
4
6
8
10
12
14
16
18
20
22
24
2
Accumulated Response of U to Shock1
4
6
8
10
12
14
16
18
20
22
24
Accumulated Response of U to Shock2
4
4
3
3
2
2
1
1
0
0
-1
-1 2
4
6
8
10
12
14
16
18
20
22
24
2
4
6
8
10
12
14
16
18
20
22
24
Como apreciamos en el resultado, este tipo de restricciones es útil para identi…car sistemas Var, cuya imposición se re‡eja en las funciones de impulsorespuesta acumuladas.
10.4. Midiendo los choques de Política Monetaria: Bernanke y Mihov (1998) (en progreso...) 10.4.1. Introducción Como se señaló en el capítulo anterior, este modelo representa el mercado interbancario donde opera el banco central. En el modelo se tienen dos grupos de variables, aquellas variables que no son de política, Y t y las variables de política, P t . El modelo puede escribirse entonces como: k
Y t =
X i=0
Bi Y ti +
k
X i=0
C i P ti + Ay V t
132CAPíTULO 10 MODELOS VAR ESTRUCTURALES EN LA PRÁCTICA k
P t =
X i=0
k
Di Y ti +
X i=0
Gi P ti + AP V tP
El conjunto de variables Y t esta dado por, el producto bruto interno, el de‡actor implícito del GDP, y un índice de precios de comodities: Y t= [GDP t ; t ; W P It ]
El conjunto de variables de política esta conformado por la tasa de interés de política, F F , el nivel de reservas no prestadas, NBR y el nivel de reservas totales j.T R P t = [FF;NBR;TR]
El nivel de reservas prestadas corresponde a las reservas generadas a través de la ventanilla de descuento del banco central, mientras que las no prestadas, a las generadas por operaciones de mercado abierto realizadas por el banco central. En forma simpli…cada, el modelo puede escribirse como: A
0@ 1A Y t
P t
= C (L)
0@ 1A 0@ 1A Y t1
P t1
+B
vtY
vtP
donde A representa las relaciones contemporáneas entre todas las variables del sistema y C (L) es la matriz de polinomios de rezagos, conformable
0@ 1A
con el vector de variables del sistema (de orden …nito). Luego, v
vtY
vtP
es el vector de choques estructurales, tanto de variables que son de política, como de aquellas que no lo son. La existencia de elementos distintos de cero fuera de la diagonal de B permite que algunos choques afecten directamente a más de una variable en el sistema. Este aspecto es clave para los propósitos de estimación del modelo, dejando claro que hace con …nalidad de efectuar un posterior análisis de politica. Modelos anteriores podian tener una especi…cacion similar, sin embargo eran elaborados asumiendo total exogeneidad por parte de las variables de política (P t ). Ello implicaba que el unico propósito del modelo era obtener multiplicadores dinámicos de política, señalando una relación de causalidad directa entre instrumentos y variables macroeconómicas. Dicha identi…cación era bastante limitada, porque como sabemos, el manejo de instrumentos de política se encuentra condicionado al comportamiento de la economía
10.4 MIDIENDO LOS CHOQUES DE POLíTICA MONETARIA: BERNANKE Y MIHOV (1 agregada y no es completamente exógeno. Esta modalidad, conocida como la comisión de Cowles, fue descartada para evaluación de medidas de política. Los modelos VAR en cambio, son estimados con el propósito de obtener evidencia empírica sobre la respuesta de variables macroeconómicas ante impulsos de politica monetaria, de tal forma que es posible discernir sobre cuál es el modelo teórico (especi…cación estructural) que mejor describe la realidad estudiada. Es por eso muy importante el conjunto de restricciones que se aplican al modelo para identi…car acciones de política, teniendo en cuenta la potencial endogeneidad de los instrumentos utilizados.
10.4.2. Identi…cación y estimación Naturalmente, sabemos que el modelo estructural no es directamente observable. Como se vo anteriormente, lo que observamos es la forma reducida. Así tenemos:
0@ 1A Y t
P t
= A 1 C (L)
0@ 1A 0@ 1A Y t1
+
P t1
uY t
uP t
donde u representa los errores de la forma reducida, vector distribuido normal e independientemente, con matriz de covarianzas . La relacion entre estos errores y los de la forma estructural es como sigue: A
0@ 1A 0@ 1A uY t
uP t
= B
vtY
vtP
que en forma simpli…cada es:
ut = A 1 Bv t
De esta última expresión se deriva la relación entre matrices de covarianzas:
0
0
0
E ut ut = A 1 BE vt vt B A1
En forma empírica tenemos:
^
^ 1 ^ ^
=A
BI B
0
^ 1
A
0
0
134CAPíTULO 10 MODELOS VAR ESTRUCTURALES EN LA PRÁCTICA
Descomposicion de Choleski La tradicional descomposición triangular será de la siguiente forma:
A =
0B BB @
1
0 0
:::
0
a21 1 0
:::
0
::: ::: :::
:::
:::
an1 ::: ::: ann1 1
1C CC A
,
B=
0B BB @
b11 0 ::: 0 0 b22 ::: 0 ::: ::: ::: ::: 0
0 ::: bnn
1C CC A
De esta forma el sistema quedaría exactamente identi…cado, donde el resultado del modelo dependerá del ordenamiento inicial de las variables (por exogeneidad). En otras palabras, se trata de una estructura recursiva.
Modelo Estructural con restricciones contemporáneas En este caso la …jación de restricciones en las matrices A y B será en base información económica sobre cómo se relacionan las variables que integran el sistema. Desde esta perspectiva, lo que se hace es incorporar elementos de teoría económica a un modelo de series de tiempo completamente ateórico y puramente matemático. La información adquirida fuera del modelo puede ser útil para que este último replique la realidad económica que estamos tratando de representar, lo cual resulta completamente válido. De esta forma, las funciones de impulso respuesta desprendidas de este modelo tendrán un comportamiento coherente, derivado de un conocimiento previo. El obejtivo en este caso es identi…car el mecanismo de transmisión de política monetaria, por lo que las restricciones de identi…cación deberán estar fundadas en ’hechos estilizados’ y no en predicciones propias del modelo. El sistema VAR cuenta con las seis variables señaladas anteriormente: El conjunto de variables Y t esta dado por, el producto bruto interno, el de‡actor implícito del GDP, y un índice de precios de comodities: Y t= [GDP t ; t ; W P It ]
El conjunto de variables de política esta conformado por la tasa de interés de política, F F , el nivel de reservas no prestadas, N BR y el nivel de reservas totales, T R. 0
P t = [FF;NBR;TR]
Los supuestos principales de identi…cación son:
10.4 MIDIENDO LOS CHOQUES DE POLíTICA MONETARIA: BERNANKE Y MIHOV (1 - Ortogonalidad en las perturbaciones estructurales. - Las variables macroeconómicas Y t , no reaccionan contemporáneamente ante cambios en las variables de política P t , pero no viceversa. - Las restricciones impuestas en el bloque monetario del modelo re‡ejan el modo de operar de la política monetaria. Las relaciones contemporáneas entre la tasa de la Fed y el mercado de reservas se derivan a continuación:
Parte II
Modelos de Series de Tiempo no Estacionarias
137
10.5 INTRODUCCIÓN
10.5. Introducción
139
Capítulo 11
Raíces Unitarias y Cointegración 11.1. Introducción El análisis de Series Económicas hoy en día es un tema de alta relevancia, en especial si se desea emitir juicios de valor relacionados con el comportamiento de las mismas. Es por esta razón que estudiamos la presencia de raíces unitarias en las series con las que trabajamos, en vista de tomar acciones correctivas en caso estas se presenten. La presencia de una raíz unitaria en una serie de tiempo económica revela directamente inestabilidad en la misma, siendo necesario el aplicar una transformación a la misma en vista de que sea posible trabajar con ella mediante los métodos estudiados anteriormente. En estricto, una serie integrada o no estacionaria es per se inestable en el tiempo, lo cual implicaría que si sufre una perturbación temporal el efecto en ella será permanente. Por tal motivo, el efectuar regresiones ó predicciones con una serie integrada puede llevarnos a emitir conclusiones erradas, situación que puede evitarse en la medida en que se lleve a cabo un proceso adecuado de identi…cación de este fenómeno. De acuerdo a la descripción anterior podemos entender que una serie económica es susceptible de ser integrada y es más, dicha característica debe tener algún signi…cado económico. En términos del análisis económico, lo que se busca es determinar si una variable reacciona de forma temporal ó permanente frente a los choques que se dan normalmente en la economía nacional e internacional. Variables, por ejemplo, como el producto bruto (en niveles) ó la Tasa de In‡ación anualizada (en porcentajes) se dice que 141
142 CAPíTULO 11 RAíCES UNITARIAS Y COINTEGRACIÓN pueden ser integradas en la medida en que son afectadas tanto por choques de oferta ó demanda y que no necesariamente pueden regresar a la tendencia que presentaban antes de la ocurrencia de dicha perturbación.
11.2. Tests de Raíces Unitarias1 Los test de raíces unitarias son utilizados para determinar si una serie integrada es estacionaria en tendencia o en diferencias. Es decir, permiten dar un aproximado del verdadero proceso generador de datos. Dadas las diferencias estadísticas entre una serie integrada y una estacionaria, usualmente se a…rma que una simple inspección visual permite identi…car el grado de la misma, así como el proceso generador de datos. Un primer ejemplo nos ilustrará acerca de este tema: Suponga dos procesos estocásticos tales como: yt =
8< :
yt1 + "t 0;9yt1 + "t
, zt =
8< :
0;05 + zt1 + "t 0;05t + "t
En ambos casos el primer proceso corresponde a un camino aleatorio, siendo incorporando una deriva en el segundo. Por su parte, el segundo proceso de (y) corresponde a un AR(1) estacionario y el segundo proceso de (z) corresponde a un proceso estacionario en tendencia. Suponga entonces que se realiza el experimento para procesos de 100 observaciones, buscando evaluar grá…camente los resultados. El código en Eviews para este experimento se presenta a continuación: ’Raíces Unitarias ’Ejemplos Chumacero (2001) ’Fernando Pérez Forero ’19/01/05 ’===================================================== !T=100 work…le u !T ’AR(1) !rho=0.9 smpl @all series e=0.01*nrnd !ci=0.0001 1
Tomado de Chumacero (2001)
11.2 TESTS DE RAíCES UNITARIAS 2 smpl 1 1 series y=!ci smpl @…rst+1 @last series y=!rho*y(-1)+e smpl @all genr ar1=y dy ’Random walk !rho=1 smpl 1 1 series y=!ci smpl @…rst+1 @last series y=!rho*y(-1)+e smpl @all genr rw=y dy group rvsar ar1 rw show rvsar.line ’trend !beta=0.05 smpl @all series e=0.01*nrnd !ci=0.000 series y=!beta*@trend+e genr trend=y dy ’Random walk + drift !alpha=0.05 !rho=1 smpl 1 1 series y=!ci smpl @…rst+1 @last series y=!alpha+!rho*y(-1)+e smpl @all genr rwd=y dy group rw_trend rwd trend show rw_trend.line
143
144 CAPíTULO 11 RAíCES UNITARIAS Y COINTEGRACIÓN .08
5
4
.04 3
.00
2
1
-.04 0
-.08
-1
10
20
30
40
50
AR1
60
70
80
90
100
10
RW
20
30
40 RWD
50
60
70
80
90
100
TREND
Observando los grá…cos, ¿es posible determinar si una serie es integrada?
11.3. El Test de Dickey – Fuller aumentado (ADF) Tradicionalmente los tests de raíces unitarias son utilizados para evaluar el verdadero proceso generador de datos (series estacionarias en media, tendencia, etc.). Como ha sido comprobado, la identi…cación grá…ca no es del todo convincente, lo cual nos lleva a buscar un soporte numérico y estadístico. El test más popular para llevar a cabo este procedimiento es el de Dickey-Fuller aumentado (ADF), que viene dado por: p
yt = + t + yt1 +
X i=1
i yti + "t
donde denota el número de rezagos necesario para que "t sea ruido blanco. La ecuación anterior es estimada por MCO con el objetivo de contrastar la Hipótesis H 0 : = 1 (yt tiene una raíz unitaria) mediante el ^ estadístico t del estimador . Tal como es documentado por Dickey y Fuller (1979) 3 , la distribución asintótica de este estimador no es estándar y por tal razón, si se utilizan los métodos tradicionales de contraste de hipótesis se producirá un error de Tipo I, Chumacero (2001) 4 detalla este problema grá…camente.Considere el siguiente proceso estocástico: 3
Dickey y Fuller (1979): Distributions of the estimates for autorregressive time series with a unit root. JASA, 74, 427-431. 4 Chumacero, R. (2001): "Testing for unit roots using economics", Documentos de Trabajo del Banco Central de Chine 102, Banco Central de Chile.
11.3 EL TEST DE DICKEY – FULLER AUMENTADO (ADF) 145 yt = 0;1 + yt1 + "t ; "t
N (0; 1)
Suponga que esta serie es generada 1000 veces bajo las mismas condiciones, donde cada muestra es también de 1000 observaciones. Nos interesa comparar la distribución del estadístico t para el test de Dickey – Fuller aumentado contra la de una normal, para lo que estimamos este valor en cada réplica del experimento. Luego gra…camos los histogramas para cada una de las distribuciones.El código en Eviews para este experimento se presenta a continuación: ’Ejemplo 2 Chumacero(2001) ’Fernando Pérez Forero ’19/01/05 ’======================================= !T=1000 work…le u !T !exp=1000 vector(!exp) td ’Random walk + drift for !i=1 to !exp smpl @all series e=nrnd !ci=0.0001 !alpha=0.05 !rho=1 smpl 1 1 series y=!ci smpl @…rst+1 @last series y=!alpha+!rho*y(-1)+e smpl @all genr rwd=y dy freeze(temp) rwd.uroot(adf,const,lag=3,save=mout) td(!i)=temp(7,4) d temp d mout next smpl @all mtos(td,tds) freeze(hist_norm) e.hist freeze(hist_tds) tds.hist
146 CAPíTULO 11 RAíCES UNITARIAS Y COINTEGRACIÓN 140
120
100
80
60
40
B
20 C A
0 -4,0
-3,4
-2,8
-2,2
-1,6
-1,0
-0,4 norm
0,1
0,7
1,3
1,9
2,5
3,1
4,0
t
Las líneas verticales corresponden a los valores críticos para cada una de las observaciones. El área denotada por “A” corresponde a la región de rechazo para un test asintóticamente normal. Asimismo, el área denotada por “C” corresponde a la región de rechazo para el test ADF. Por tanto, el área “A+B” representa la distorsión que genera el utilizar un valor crítico estándar para este tipo de test. Dado este problema, es bastante frecuente declarar una serie como estacionaria cuando en realidad es integrada. En particular, es importante señalar que el ejemplo anterior se hizo con la …nalidad de ilustrar los posibles problemas que se pueden encontrar al utilizar el test ADF, tales como utilizar los valores críticos incorrectos, así como identi…car una serie como integrada cuando en realidad es estacionaria, pero cercana a la raíz unitaria. No obstante, el test descrito es llevado a la práctica en el apartado siguiente.
11.4. El Test de Dickey – Fuller aumentado (ADF) en Eviews En Eviews lo importante es indicar el proceso generador de datos antes de ejecutar el test. Son tres ecuaciones las que se estiman alternativamente:
11.4 EL TEST DE DICKEY – FULLER AUMENTADO (ADF) EN EVIEWS 147 Estacionariedad
Test ADF p
Estacionaria en tendencia
X X X
yt = + t + yt1 +
p
Estacionaria en media ( = 0)
yt = + yt1 +
p
Autorregresivo ( = = 0)
yt = yt1 +
i=1
i=1
i=1
i yti + "t
i yti + "t
i yti + "t
El orden correcto es el de estimar la ecuación estacionaria en tendencia ^ y contrastar el valor de . Si en caso el resultado es una serie estacionaria ( < 1), el proceso termina. Caso contrario, el siguiente paso es el de veri…car la signi…cancia de . En caso no sea signi…cativo, deberá evaluarse la ^ ecuación estacionaria en media y volver a contrastar el valor de . Si en caso el resultado es una serie estacionaria, el proceso termina. Caso contrario, el siguiente paso es el de veri…car la signi…cancia de . En caso no sea signi…cativo, deberá evaluarse la ecuación autorregresiva y volver a contrastar el ^ valor de . En caso el resultado es una serie estacionaria, el proceso termina, caso contrario la serie tiene una raíz unitaria. El código en Eviews para este experimento se presenta a continuación: ’ Test ADF Secuencial ’ Fernando Pérez Forero ’ 01/02/06 ’ ===================================================== ’ Argumentos: ’ Serie: %0 ’===================================================== ’ Resultados (ver trend, inter, nin): ’ 0: PGD Inválido ’ 1: Serie Integrada y PGD válido ’ 2: Serie Estacionaria y PGD válido ’===================================================== genr y= {%0} ’Tendencia freeze(trend_adf) y.uroot(adf,t) %data=trend_adf(3,1) %lag=@mid( %data,13,2) !lag=@val( %lag) !tstat=@val(trend_adf(7,4))
148 CAPíTULO 11 RAíCES UNITARIAS Y COINTEGRACIÓN !MKt=@val(trend_adf(9,4)) !Bt=@val(trend_adf(24+!lag+2,5)) vector(1) tend if abs(!tstat)
0.05 then tend(1)=0 ’Intercepto freeze(inter_adf) y.uroot(adf,c) !tstat=@val(inter_adf(7,4)) !MKt=@val(inter_adf(9,4)) !Ct=@val(inter_adf(24+!lag+1,5)) vector(1) inter if abs(!tstat)0.05 then inter(1)=0 ’Ninguno freeze(nin_adf) y.uroot(adf,n) !tstat=@val(nin_adf(7,4)) !MKt=@val(nin_adf(9,4)) vector(1) nin if abs(!tstat)
Replicaremos este procedimiento para la tasa de crecimiento del PBI Peruano, tomando datos mensuales y desde enero de 2002. Intuitivamente esta serie debería ser estacionaria, debido a que los choques en la economía
11.4 EL TEST DE DICKEY – FULLER AUMENTADO (ADF) EN EVIEWS 149 agregada y los quiebres que esta implica se re‡ejan en el nivel de producto (boom de producción, crisis …nanciera internacional, etc.), pero en el largo plazo estas variaciones son relativamente constantes. Grá…camente tenemos:
10
8
6
4
2
0 2002
2003
2004
2005
VPBI
Una simple inspección nos lleva a a…rmar que esta serie es estacionaria, con una ligera inclinación positiva, lo cual implicaría la posible presencia de un componente tendencial. De acuerdo al programa descrito anteriormente, los resultados indican que esta serie es estacionaria y que efectivamente, cuenta con un componente tendencial (signi…cativo al 10 %). Adicionalmente, no fue necesario utilizar primeras diferencias rezagadas en la estimación.
150 CAPíTULO 11 RAíCES UNITARIAS Y COINTEGRACIÓN
11.5. Raíces unitarias en presencia de cambios estructurales En términos estadísticos el aceptar una hipótesis falsa es considerado como un error de Tipo II, es decir, si concluimos que una serie es integrada cuando no lo es nos encontraremos en el escenario descrito. Cuando las series de tiempo presentan cambios estructurales y son sometidas al test ADF el resultado es susceptible de ser fallido. En otras palabras, se declarará una serie como integrada cuando no necesariamente lo es. En vista de ello, el Test de Zivot y Andrews se encarga de solucionar este problema, encontrando el verdadero proceso generador de datos y calculando el punto aproximado donde se genera el quiebre en la serie. En estricto, el test de Zivot y Andrews calcula el valor del test ADF y sus variantes para cada observación en el
11.5 RAíCES UNITARIAS EN PRESENCIA DE CAMBIOS ESTRUCTURALES 151 tiempo, colocando los resultados en un grá…co. De esta forma podremos visualizar si la serie presenta un cambio drástico en algún punto del tiempo para este estadístico.
Capítulo 12
Cadenas de Markov y Modelos de Cambios de Régimen 12.1. Introducción y motivación Durante los últimos años se ha venido haciendo énfasis en el uso de técnicas de modelación avanzadas dentro del marco del análisis de series de tiempo económicas, hecho motivado principalmente por las de…ciencias del análisis de regresión clásico en materia de robustez y poder predictivo. En efecto, en el trabajo diario temas como cambios estructurales, inestabilidad de parámetros, no estacionariedad, volatilidad, outliers, etcétera, nos llevan a cuestionar las conclusiones a las que nos hacía llegar el análisis de modelos ARIMA ó ARMA. Ante esta problemática han surgido modelos que, apuntando a obtener una transformación plausible de la serie a analizar, lograban aproximarse a la metodología en cuestión. Por otro lado, también surgieron modelos que, en vez de transformar la serie, apuntaban mas bien a extraer ó ‘…ltrar’ un componente no observable de la serie. Dicho componente es el conocido ‘estado’ de la variable y es el que guía las ‡uctuaciones principales de la serie bajo análisis. De hecho, han surgido diversos algoritmos que se encargan de estimar dicho componente en base a la información disponible y el más conocido es el denominado Filtro de Kalman. El estado estimado se desenvuelve dentro de un espacio determinado ó domino, y es por esta razón que se conoce a esta clase de modelos como de espacio de estados (state-space). Los modelos de Markov-Switching forman parte de este grupo de algo153
154CAPíTULO 12 CADENAS DE MARKOV Y MODELOS DE CAMBIOS DE RÉG ritmos que se encargan de …ltrar el componente no observable de una serie de tiempo. La diferencia radica en que en éstos modelos el estado es una variable aleatoria discreta. Es decir, el espacio del estado es un conjunto de puntos discretos asociados a un proceso estocástico determinado. En un modelo de Markov-Switching este proceso estocástico sigue una Cadena de Markov. En estricto, una Cadena de Markov es un proceso estocástico que cumple con la propiedad de Markov, que señala que la probabilidad de encontrarse en un estado determinado j en un instante t es condicionada únicamente al estado en se que ubicó el proceso en t-1. En otras palabras, la única información relevante para poder predecir donde se encontrará dicha variable en el instante t es la posición que ocupó en t-1 y mas bien el resto de información hacia atrás ya no es de utilidad. En este capítulo veremos al detalle en qué consiste una Cadena de Markov, cómo se caracteriza, cuáles son sus propiedades principales y cómo se construye dicho proceso, para luego pasar a aplicar dicho concepto en el análisis de series de tiempo mediante la simulación de un modelo sencillo R. utilizando MATLAB
12.2. Cadenas de Markov 12.2.1. De…niciones principales Proceso Estocástico: Secuencia o vector de variables aleatorias. fyt g1 t=0 Propiedad de Markov: Se dice que un proceso estocástico fxt g tiene la propiedad de Markov si 8k 2;:::;t se cumple que Pr (xt+1 jxt ; xt1 ;:::;xtk ) = Pr (xt+1 jxt ), es decir, los rezagos no añaden información a la predicción. Este proceso es caracterizado como una Cadena de Markov. Por lo tanto, una Cadena de Markov invariante en el tiempo es determinada por tres objetos: Un vector x 2 Rn que representa los valores posibles del estado de sistema. En otras palabras, el espacio del estado es de orden …nito (n < 1). Una matriz de transición P (nxn) que mide la probabilidad de moverse de un estado a otro. El elemento pij mide la probabilidad de pasar del estado i en el instante t al estado j al período siguiente (t + 1). pij = Pr (xt+1 = x j jxt = x i ). P (nxn) es entonces una matriz estocástica.
155
12.2 CADENAS DE MARKOV
Un vector 0 2 Rn que representa la distribución inicial de los estados. El elemento 0i mide la probabilidad de encontrarse en el estado i en el instante t = 0. Bajo estas propiedades podemos asumir que: n
8i = 1;:::;n se cumple que
X
pij = 1. La suma de todas las …las de
j=1
P es igual a 1, lo cual tiene sentido desde que cada …la representa el universo de estados a terminar partiendo del estado i . Nota: F il = j y C ol = i . n
X
0i = 1: El universo de estados en los cuales se puede empezar la
i=1
Cadena de Markov. La matriz estocástica P de…ne la probabilidad de moverse del estado i al estado j de un período al siguiente. Luego, si quisiéramos medir la probabilidad de pasar del estado i al estado j luego de dos períodos tendríamos lo siguiente:
Pr (xt+2 = x j xt = x i ) n
=
X h=1
j
Pr (xt+2 = x j xt+1 = x h ) P r (xt+1 = x h xt = x i )
j
j
Es decir, se tiene en cuenta el universo de puntos intermedios sobre los cuales se puede situar x en t + 1: n
Pr (xt+2 = x j xt = x i ) =
j
X
pih phj = p2ij
h=1
En general tenemos que: Pr (xt+k = x j jxt = xi ) = p kij . Probabilidad de pasar del estado i al estado j luego de k períodos.Las probabilidades vistas hasta ahora se re…eren a distribuciones condicionales, es decir, que dependen de la posición o estado anterior. Si por el contrario quisiéramos determinar la probabilidad de encontrarse en un estado j en un instante determinado y de forma incondicional debemos trabajar con lo siguiente:
156CAPíTULO 12 CADENAS DE MARKOV Y MODELOS DE CAMBIOS DE RÉG Partiendo de la distribución incondicional inicial 0 es posible determinar las distribuciones iniciales para cada período, trabajando con la matriz de transición P . 0
1 = Pr (x1 ) = 0 P
0
0
2 = Pr (x2 ) = 1 P = 0 P 2 ::: k = Pr (xk ) = k1 P = 0 P k 0
0
12.2.2. Distribuciones estacionarias En base a la ecuación anterior podemos inferir cómo hallar una distribución estacionaria incondicional, siguiendo la ley de movimiento: t+1 = t P . Luego, si la distribución es estacionaria, entonces t+1 = t = , lo cual nos llevaría a: = P =) (I P ) = 0. Tomando la transpuesta: I P = 0. Cuando hablamos de estacionariedad nos referimos a momentos invariantes en el tiempo. Se trata en este caso de estacionariedad en sentido fuerte, es decir, que todos los momentos de esa distribución son constantes en el tiempo. En términos del álgebra lineal esta ecuación de…ne un vector propio asociado a un valor propio de P igual a la unidad ( = 1). En otras palabras, para que exista una distribución estacionaria es necesario que P tenga un valor propio unitario. Luego, analizando la de…nición de P sabemos que 0
0
0
0
0
0
0
0
n
X
pij = 1, lo cual garantiza que existe al menos un ( = 1) y por tanto
j=1
existe una que satisface esta ecuación. Dependiendo de los elementos de P es posible que exista más de un ( = 1), lo cual indicaría que no es única.
12.2.3. Derivación de la probabilidades del estado estacionario Dada la matriz de transición de probabilidades P de orden M y un vector columna iM = 1M , donde todas sus componentes son iguales a 1. Entonces n
0
0
se cumple que iM P = i M desde que
X
pij = 1.Sea luego un vector M x1 de
j=1
0
probabilidades de estado estacionario, tenemos que iM = 1. De acuerdo a la
157
12.2 CADENAS DE MARKOV
24 35
de…nición anterior se cumple que I M
24
3 5 0
I M
P
0
iM
24 24
=
35 35
0M
.Sea A =
1
0
P = 0M . . Tenemos entonces: 0
I M
P
0
una matrix (M + 1) xM .
iM
Dado que A no es cuadrada, es imposible despejar a menos que tomemos el producto interno y pre-multipliquemos por la inversa, tal como sigue: 0
= A A
1 0 0M A . 1
0
1
A0 , será igual a la última columna de Z . Sea entonces Z = A A R que genera este proceso se llama markov_ss.m. La rutina de Matlab
12.2.4. Estacionariedad Asintótica Es difícil hablar de distribuciones con momentos exactamente iguales a lo largo del tiempo, por lo que empíricamente es recomendable trabajar con una suerte de estacionariedad asintótica ó en el límite. El experimento es como sigue: Partiendo de una distribución arbitraria 0 , ¿las distribuciones condicionales t convergen a una estacionaria? La pregunta puede ser respondida probando la existencia del límite: L{m t = 1 . Luego, si es que este existe es necesario también determinar t!1 si está asociado a la condición inicial 0 . Esto se puede probar realizando simulaciones sobre un mismo modelo, analizando la sensibilidad del límite a la condición inicial. Si la respuesta es negativa entonces podemos a…rmar que el proceso es asintóticamente estacionario y con una distribución única e invariante. De…nición: Sea 1 el único vector que satisface I P = 0. Si 8 0
0
t
0
se cumple que P converge siempre a 1 , decimos entonces que la Cadena de Markov es asintóticamente estacionaria con una distribución única e invariante.
12.2.5. Simulación de una Cadena de Markov Dada una matriz de transición P asociada a un proceso de markov S t , sabemos que es posible calcular la distribución estacionaria de este proceso, dada por . Si tomamos esta distribución como punto de partida 0 , procederemos a simular luego los valores de S t .
158CAPíTULO 12 CADENAS DE MARKOV Y MODELOS DE CAMBIOS DE RÉG La lógica de la simulación va de la mano con el tema de generación de números aleatorios. Para iniciar la cadena de markov es necesario tomar valores aleatorios iniciales, los cuales nos van a permitir determinar en qué estado se ubicará la cadena en base a la distribución asignada. Un ejemplo facilitará el entendimiento del proceso. Seguiremos la rutina dada en ejemplo_markov.m: Supongamos que el espacio de estados viene dado por xt = f0; 1g y que la matriz de transición asociada es P =
24 35 0;7 0;3
. Con estos datos la
0;2 0;8
distribución estacionaria es fácilmente obtenida, .
Para el tema de la generación de número aleatorios se utilizará la distribución uniforme entre 0 y 1, que es consistente con la necesidad de hallar valores en este rango, desde que trabajamos con probabilidades. ¿Cómo se determina la posición inicial en la cadena? Sea Generamos un numero aleatorio ~. Luego, si , entonces la posición queda descartada como inicial. Si mas bien , entonces es la posición inicial. En nuestro ejemplo: Si , entonces queda descartado y mas bien queda como posición inicial. El vector de estado en sería De no cumplirse la premisa, la solución es contraria a la anterior. Sea , donde es la columna i-ésima de . Para hallar la distribución de probabilidad en debemos pre-multiplicar por . De esta forma, se tiene la distribución de probabilidad de condicionada a ( ). Del mismo modo, generamos y ~encontramos la posición discriminando por el número aleatorio. La simulación se desarrolla de la misma forma durante T iteraciones (asignadas arbitrariamente). Lo más relevante es capturar la posición en cada periodo, para luego gra…carla. La rutina ejemplo_markov.m lleva a cabo este proceso, donde la simulación es hecha para T = 50, resultando el siguiente grá…co:
12.3 MODELOS DE MARKOV SWITCHING
159
12.3. Modelos de Markov Switching 12.3.1. Introducción El objetivo de este apartado es estimar modelos univariantes con cambios de régimen, donde la variable de estado del modelo es discreta y sigue una Cadena de Markov. En líneas generales el uso de las Cadenas de Markov en el análisis de series de tiempo es importante porque: Permite estimar modelos con quiebres estructurales sin conocer expresamente en qué puntos del tiempo se encuentran estos. Los cambios de régimen son endógenos a las características del modelo y variables explicativas.
160CAPíTULO 12 CADENAS DE MARKOV Y MODELOS DE CAMBIOS DE RÉG Brinda una comprensión mayor sobre la noción de ‘estado’ y lo que representa esta variable. Los primeros en introducir los modelos de Markov Switching fueron Goldfeld y Quandt (1973), quienes lo aplicaron a un modelo de varias variables explicativas. Posteriormente Hamilton (1989) extendería este tema al análisis uni y multivariado de procesos autorregresivos temporales. Desarrollaremos en este apartado un modelo AR(1) con cambios estructurales en el parámetro de persistencia, que se deriva a continuación: yt = S t yt1 + "t ; t = 1; 2;:::;T; donde "t
N (0; 2S ) t
S t) + 1S t ) Si S t = 1 ) 1 Si S t = 0 ) 0 2 2 2 ) Si S t = 1 ) 21 S = 0 (1 S t ) + 1 S t Si S t = 0 ) 20 S t = f0; 1g, Régimen 0 o 1. S t = 0 (1 t
Nótese que si S t es conocido a priori, entonces hace las veces de una variable dummy tradicional. Bajo este contexto, la función de Log-Verosimilitud viene dada por : T
Ln (L) =
X
Ln (f (yt S t )) ; dondef (yt S t ) =
t=2
j
j
El problema de maximización sería:
q
1 2 2S t
Max
exp
! yt
S t yt1 2 2S t
Ln (L)
f ; ; ; g 0
1
2 0
2 1
Las di…cultades surgen cuando el vector S t no es observado ex – ante, lo que hace que los algoritmos predictores-correctores sean los principales candidatos a ser utilizados para …ltrar dicho vector. Esta clase de algoritmos predice, con la información disponible en t 1 , el valor de S t con cierta con…abilidad y luego, al realizarse yt , esta información es actualizada con una ganancia.
12.3.2. Estimación del modelo Paso 1: Considere la densidad conjunta de y t y S t (no observado), que puede ser expresada como el producto de las densidades marginales: f yt ; S t jt1 = f yt jS t ; t1 f S t jt1 , donde t1 representa el set de información disponible en el instante t 1, que para …nes prácticos resulta siendo el conjunto de rezagos de y t .
2
161
12.3 MODELOS DE MARKOV SWITCHING
Paso 2: Para obtener la densidad marginal de yt debe integrarse con respecto a S t , es decir, sumar las densidades de yt para todos los posibles valores de S t . 1
j X j Xj j j p p X ( X j j )
f yt t1 =
f yt ; S t t1
S t =0 1
=
=
f yt S t ; t1 f S t t1
S t =0 1 2 20
exp
fyt
y 2
0 t1 2 0
g
2
Pr S t = 0 t1 +
Por lo tanto, la función de Log-verosimilitud viene dada por: T
Ln (L) =
exp
1
Ln
t=2
1 2 21
,
f yt S t ; t1 f S t t1
S t =0
que puede ser interpretada como un promedio ponderado de densidades marginales. Entonces, el problema está en calcular correctamente las ponderaciones ó probabilidades de S t . Esto puede ser realizado fácilmente mediante un …ltro que presentaremos a continuación. Sin embargo, como toda estimación no lineal, es necesario establecer una condición inicial para que esta pueda llevarse a cabo. En otras palabras, debemos suponer un comportamiento a priori de S t plausible. Por simplicidad, asumiremos que S t sigue una Cadena de Markov de primer orden, lo cual implica que solamente el primer rezago es relevante para predecir la posición futura de la variable. Es importante señalar que los procesos de Markov no son la única alternativa para construir modelos de saltos discretos (switching), pero en este caso consideramos que es lo que más se aproxima al análisis que se desea realizar. Construimos entonces las probabilidades de la función de verosimilitud: Pr (S t = 0 S t1 = 0) = p = Pr (S t = 1 S t1 = 1) = q =
j j
1 1+exp( p0 ) 1 1+exp(q0 )
con sus respectivos complementos. Estas probabilidades harán las veces de las ponderaciones de las densidades marginales de yt en la función de verosimilitud. Fijando valores iniciales para = 0 ; 1 ; 20 ; 21 ; p0 ; q 0 es posible comenzar a iterar la función de verosimilitud a …n de obtener los valores óptimos R se llama asociados a los datos de yt . La función construida en Matlab Loglike_ar1.m y presenta justamente estos seis argumentos. El problema de maximización sería ahora: MaxLn (L). La rutina de
R que ejecuta este proceso es ar1.m. En estricto, el programa minMatlab
fyt
y 2
1 t1 2 1
2
g Pr S = 1j t t 1
162CAPíTULO 12 CADENAS DE MARKOV Y MODELOS DE CAMBIOS DE RÉG imiza el negativo de esta función, que para …nes prácticos nos brinda el mismo resultado.
12.3.3. Filtrado de probabilidades ^ ^ ^ 2 ^ 2 ^ ^
Una vez estimados 0 ; 1 ; 0; 1 ; p; q es posible …ltrar las probabilidades de S t condicionales a la información disponible en t 1. Dado que es una variable no observada, adoptamos el siguiente algoritmo para el cálculo de las probabilidades.
Paso 1: Dada Pr S t1 = i j t1 , i = 0; 1 al inicio de la iteración t , las ponderaciones Pr S t = j j t1 se calculan como sigue:
1
X X
Pr S t = j t1 =
j
i=0
Pr S t = j; S t1 = i t1
j
1
=
i=0
Pr (S t = j; S t1 = i) Pr S t1 = i t1
j
Observamos que el primer término corresponde a la probabilidad de transición. Asimismo, el segundo es la probabilidad observada al inicio de la iteracion, por lo que ambos elementos son conocidos. A este proceso se le conoce como predicción de probabilidades. Paso 2: Una vez que llega la realización de yt (al término de la iteración t), podemos actualizar la probabilidades como sigue:
X
Pr (S t = j t ) = Pr S t = j yt ; t1 =
j
=
j
f (yt S t = j;t
1
j
i=0
1
j j
1
f (yt t
1 ) Pr(S t = j j t 1 )
f (yt S t =i;t
j
f (S t = j yt ;t 1
)
)
) Pr(S t =ijt 1 )
donde t = t1 ; yt . Esta probabilidad actualizada sirve como dato para la iteración siguiente. Se asume que la matriz de transición es invariante en el tiempo. Es necesario iniciar el …ltro de algún modo, para lo cual lo más recomendable sería tomar la distribución estacionaria en base a la matriz de transición R markov_ss.m estimada en el apartado anterior. El programa de Matlab realiza este proceso, tal como fue indicado líneas arriba. Del mismo modo, la función Filter.m realiza el cálculo de la probabilidades tanto predichas como actualizadas. La rutina ar1.m llama también a esta función luego de obtener los parámetros estimados en el apartado anterior.
163
12.3 MODELOS DE MARKOV SWITCHING
12.3.4. Suavizado (Smoothing) de probabilidades Las probabilidades …ltradas en el paso anterior eran condicionadas a la información disponible en t 1 (predicha) y en t (actualizada). El algoritmo propuesto por Kim se encarga de calcular las probabilidades en cada instante t condicionadas a la información disponible de toda la muestra ( T ). Dicho de otro modo, se ajustan las probabilidades dado que existe una mayor información disponible, lo que para …nes prácticos implica que es posible conocer la probabilidad de encontrarse en un estado j en el instante t , dado que el estado se ubicó en k en el instante t + 1. La probabilidad conjunta vendría dada por: Pr (S t = j; S t+1 = k T ) = Pr (S t+1 = k T ) P r (S t = j S t+1 = k; T )
j
j j Pr(S = j;S =kj ) = Pr (S t+1 = k jT ) Pr(S =kj ) =k;S = j)Pr(S = j j ) = Pr (S t+1 = k jT ) Pr(S Pr(S =kj ) t
t+1
t
t+1
t
t
t+1
t
t+1
t
t
Nótese que el segundo término del lado derecho viene formado por las probabilidades …ltradas antes y también por las probabilidades de transición estimadas. 1
Luego: Pr (S t = j j T ) =
X
Pr (S t = j; S t+1 = k T )
j
k=0
Implementando el …ltro: Al alcanzar la iteración T en el proceso de …ltrado se obtiene la siguiente información: Pr (S T = 0j T ), Pr (S T = 1j T ) (actualizadas) y Pr S T = 0jT 1 , Pr S T = 1jT 1 (predichas). En T 1:
Pr (S T 1 = j; S T = k T ) = Pr (S T = k T )
j
1
Pr (S T 1 = j T ) =
j
X k=0
j
Pr(S T =k S T
j
1
= j) Pr(S T
1
Pr(S T =k T
j
1
)
= j T
j
1
Pr (S T 1 = j; S T = k T )
j
R realiza este procedimiento. Los arLa función Smooth.m de Matlab gumentos en este caso son las probabilidades predichas y actualizadas calculadas por Filter.m, así como el vector de parámetros estimados por Loglike_ar1.m. La rutina ar1.m llama a esta función luego de ejecutar el …ltrado de las probabilidades. En línea con lo explicado en este apartado, podemos realizar un experimento para 15 observaciones con ar1.m, obteniendo el siguiente resultado para = 0 ; 1 ; 20 ; 21 ; p ; q :
)
164CAPíTULO 12 CADENAS DE MARKOV Y MODELOS DE CAMBIOS DE RÉG
165
12.4 ALGUNAS EXTENSIONES
12.4. Algunas extensiones 12.4.1. Aproximación discreta Generalmente se tienen muchas di…cultades a la hora de resolver modelos dinámicos, puesto que no siempre es posible obtener soluciones analíticas. Alternativamente existe un procedimiento que permite aproximar de forma discreta un proceso autorregresivo, presentando un espacio de estados y una matriz de transición. Para …nes prácticos, decimos que la solución aproximada sigue una Cadena de Markov invariante de primer orden. Tauchen (1985) desarrolló este tema tanto para modelos univariados como multivariados, dejándonos la enseñanza de que debemos aprovechar las ventajas computacionales de nuestros tiempos en vista de llegar a soluciones rápidamente y con sentido. En estricto, Tauchen …ltra los posibles valores de la solución y los agrupa en un espacio de estados, brindando también la matriz de transición asociada. En esta oportunidad solamente veremos el caso escalar. Supongamos que se tiene un proceso AR(1) tal que: yt = yt1 + " t , donde "t es un ruido blanco con varianza 2" . Sea F la función de densidad acumulada de "t tal que: Pr ("t u) = F u" . Por simplicidad se asume " = 1. Luego sea yt el proceso discreto que aproxima a y t y sean y 1 < y2 < ::: < yN los N valores equidistantes que yt puede tomar. Una alternativa para seleccionar estos valores es tomar a y N como un múltiplo m de la desviación
e
e
2" 1 2
1 2
estándar incondicional del proceso y = . Luego y1 = yN y el resto de valores son equidistantes. Con el espacio de estados aproximado, se procede a estimar la matriz de transición pij = Pr yt = y j jyt1 = y i que sigue el proceso. Sea el paso w = y j y j 1 , luego para cada i, si 2 j N 1, entonces:
e e !
pij = Pr y j
w 2
yi + "t
y j +
w 2
despejando en función a " t tenemos por diferencia de áreas: pij = F
y j +
w 2
"
y i
F
y j
w 2
y i
"
:
Caso contrario en tendremos: pi1 = F
y1 +
w 2
"
y i
opiN = 1
F
yN
w 2
"
y i
166CAPíTULO 12 CADENAS DE MARKOV Y MODELOS DE CAMBIOS DE RÉG La racionalidad para asignar las probabilidades de esta forma radica en que yt es el proceso discreto que aproxima a yt , por lo tanto, de acuerdo a la ley de movimiento de y t ) y j yi + "t . Dado que y i es …jo, entonces la distribución de probabilidad de y j es una función lineal de la distribución de "t . R tauchen.m describe este proceso en forma detalLa rutina de Matlab lada.
e
12.4.2. Modelado de preferencias Sea un consumidor con preferencias E
(X ) 1
t
t=0
C t1 1
, donde 0 < < 1
y > 1. Sea que el consumo sigue una Cadena de Markov de primer orden, el espacio de estados que puede tomar el consumo es: C t = fC 1 ; C 2 ;:::;C N g . Denotaremos P como la matriz de transición asociada a este proceso y 0 como la distribución inicial. 1
Sean ui = C 1i , vi = E escalar V = E fvi g. Sea v =
1
X
t
P t u =
t=0
1
X
(X ) (X 1
t
1
C t 1
t=0
E
1
t=0
t
1
C t 1
)
jC 0 = C i y el
(P )t u.
t=0
Luego, es fácil demostrar que: v = (I N P )1 u. Mientras no se conozca el estado inicial, el valor esperado de v viene dado por: N
V = E v =
fg
X
vi 0i .
i=1
Supongamos que el consumo tiene dos estados C = (1; 5) y que la probabilidad de comenzar en un estado o en otro está dividida en partes iguales, es decir 0 = (0;5; 0;5). Consideremos ahora dos procesos de markov con matrices de transición: P 1 =
0@ 1A 1 0 0 1
y P 2 =
0@ 1A 0;5 0;5
, manteniendo constante tanto el
0;5 0;5
espacio de estados como la distribución inicial. Sean = 0;95 y = 2;5. La idea sería determinar qué proceso de markov nos lleva a una función de valor más alta, para luego establecer un ordenamiento de preferencias entre matrices de transición.
167
12.5 REFERENCIAS
R consumo.m lleva a cabo este proceso y lo intereLa rutina de Matlab sante de este resultado es que V P = V P = 7;2630, lo que indica que el consumidor es indiferente a las matrices de transición dadas. Cabe señalar que este resultado se mantiene inclusive si llegase a variar. La interpretación de este resultado es que, dado que el consumidor tiene completa incertidumbre sobre en qué estado se va a situar inicialmente, la distribución inicial asigna 50 % para cada posibilidad. En ese instante, sobre P 1 el consumidor tiene 50 % de probabilidad de tener un consumo de 5 toda tu vida y 50 % de tener un consumo de 1, por lo que está completamente indiferente sobre qué estado alcanzará. Del mismo modo, bajo P 2 el consumidor es indiferente al inicio y a lo largo de toda su vida sobre qué estado alcanzará. De acuerdo a este planteamiento, sea cualesquiera de las dos matrices, la racionalidad del consumidor es la misma, es decir, incertidumbre total sobre qué nivel de consumo alcanzará, lo que lleva a una función de utilidad esperada igual en ambos casos. Un análisis mucho más profundo nos dice que si hacemos variar la distribución de probabilidad inicial, el consumidor estará revelando que tiene cierta información que lo hacen pensar que alcanzará un nivel de consumo determinado, en vez del otro. Si esto ocurriese entonces el consumidor podrá entonces discriminar entre un proceso de markov y otro. En breve, si el consumidor …ja una probabilidad más alta para 1, entonces preferirá P 2. Caso contrario, preferirá P 1 . Es lógico desde que si sabe que es más probable que su consumo sea bajo, no va a preferir una matriz de transición donde no va a poder salir de ese estado. Si por el contrario, sabe que es más probable que su consumo sea alto, va a preferir entonces una matriz de transición donde se asegure de no salir de ese estado. Este ejemplo es importante porque permite …jar una línea de base para la construcción de modelos de equilibrio general con cambios de régimen. 1
2
12.5. Referencias Aoki, M.(?): “Not More So: Some Concepts Outside the DSGE Framework”. Department of Economics, University of California, Los Angeles. Fischer, A.; Dueker, M & Dittmar, R. (2002): “Stochastic Capital Depreciation and the Comovement of Hours and Productivity”. Swiss National Bank and CEPR.
168CAPíTULO 12 CADENAS DE MARKOV Y MODELOS DE CAMBIOS DE RÉG Hamilton, J. (1989): “A New Approach to the Economic Analysis of Nonstationary Time Series and Business Cycle.” Econometrica, Vol.57, No 2 (Mar. 1989), 357-384. Humala, A. (2004): “Modelling interest rate pass-through with endogenous switching regimes in Argentina”. University of Warwick, June 2004. Kim, C. & Nelson, C. (1999): “State-Space Models with Regime Switching”. Capítulo 4. The MIT Press. Krolzig, H.-M. (1997): “Markov-Switching Vector Autoregressions. Modelling, Statistical Inference and Application to Business Cycle Analysis”, Lecture Notes in Economics and Mathematical Systems, Volume 454, Berlin: Springer. Ljungqvist, L & Sargent, T. (2000): “Recursive Macroeconomic Theory”. Capítulo 1. The MIT Press. Sims, C. & Zha, T. (2004): “Were There Regime Switches in U.S. Monetary Policy?” Working paper 2004-14, June 2004. Federal Reserve Bank of Atlanta. Tauchen, G. (1985): “Finite State Markov-Chain approximations to univariate and Vector Autorregressions”. Economics Letters 20 (1986) 177-181.
12.6. Apéndice: Programas de Matlab 12.6.1. ar1.m % Procedimiento que replica a Hamilton (1989) % Basado en Kim & Nelson (1998) % Markov Switching AR(1) % % Fernando Pérez Forero % 02/10/05 % clear P=[0.80,0.20;0.30,0.70]; phi1=log((1/P(1,1))-1);
12.6 APÉNDICE: PROGRAMAS DE MATLAB phi2=log((1/P(2,2))-1); beta0=0.99; beta1=0.65; sigma0=1; sigma1=0.5; T=15; rho=0.8; mu=0; y(1)=randn(1); for i=2:T y(i)=mu+rho*y(i-1)+randn(1); end y=y’; global y; theta(1)=beta0; theta(2)=beta1; theta(3)=sigma0; theta(4)=sigma1; theta(5)=phi1; theta(6)=phi2; Loglike_ar1(theta); % Seleccionando el vector de opciones default options = optimset(’fminsearch’); % Ajustando las opciones options.MaxFunEvals= 30000; options.LargeScale=’o¤’; options.MaxIter= 15000; [theta_…n,LogF,exit‡ag,output] =fminsearch(’Loglike_ar1’,theta,options); theta_…n; for i=1:4 Result(i)=theta_…n(i); end P1est=1/(1+exp(theta_…n(5))); Q1est=1/(1+exp(theta_…n(6))); Result(5)=P1est; Result(6)=Q1est; Result
169
170CAPíTULO 12 CADENAS DE MARKOV Y MODELOS DE CAMBIOS DE RÉG Pest=[P1est,(1-P1est);(1-Q1est),Q1est] LogF=-LogF exit‡ag %====================================================== % Filtering [Filt_pr,Filt_act]=Filter(theta_…n); …gure(1) subplot(2,2,1);plot(Filt_pr(:,1),’b’);title(’Pr(St=0 jt-1)’);xlabel(’t’);ylabel(’Pr(St=0)’); subplot(2,2,2);plot(Filt_pr(:,2),’r’);title(’Pr(St=1 jt-1)’);xlabel(’t’);ylabel(’Pr(St=1)’); %====================================================== % Smoothing % Tenemos Pr(ST=0jphi(T)) y Pr(ST=1jphi(T)) % Iterando hacia atrás: Smooth_m=Smooth(Filt_pr,Filt_act,theta_…n); subplot(2,2,3);plot(Smooth_m(:,1),’b’);title(’Pr(St=0 jT)’);xlabel(’t’);ylabel(’Pr(St=0)’); subplot(2,2,4);plot(Smooth_m(:,2),’r’);title(’Pr(St=1 jT)’);xlabel(’t’);ylabel(’Pr(St=1)’); %====================================================== % Simulación de Cadenas de Markov % Simulación sobre la muestra % En este caso la distribución inicial será la del estado estacionario. pi0=markov_ss(Pest); pi0=pi0’; x=[0;1]; [simul,estado1]=simulate_markov(x,Pest,pi0,T); …gure(2) subplot(2,2,1:2);plot(simul,’bo-’);title(’Cadena Simulada sobre la muestra’);xlabel(’t’);ylabel(’S(t)’); % Simulación hacia adelante % En este caso la distribución inicial será formada por las probabilidades % …ltradas para t=T (actualizadas). pi0_s=[Smooth_m(T,1),Smooth_m(T,2)]; h=20; [simul,estado2]=simulate_markov(x,Pest,pi0_s,h); subplot(2,2,3:4);plot(simul,’bo-’);title(’Cadena Simulada h períodos hacia adelante’);xlabel(’h’);ylabel(’S(T+h)’);
12.6.2. Loglike_ar1.m function LogL=Loglike_ar1(theta)
12.6 APÉNDICE: PROGRAMAS DE MATLAB
171
% Función de Log-verosimilitud para un AR(1) de T observaciones con 2 % estados. Se asume que el modelo no tiene intercepto. % theta: Vector de parámetros. % Fernando Pérez Forero % 02/10/05 global y; y0=randn(1); p1=1/(1+exp(theta(5))); q1=1/(1+exp(theta(6))); P=[p1,(1-p1);(1-q1),q1]; Pi=markov_ss(P); T=length(y); %t=1 LogL=0; alpha = p1*Pi(1,1)+(1-p1)*Pi(2,1); beta = (1-q1)*Pi(1,1)+q1*Pi(2,1); phi=normpdf(y(1),y0*theta(1),theta(3))*alpha+normpdf(y(1),y0*theta(2),theta(4))*beta; alpha_pr=(normpdf(y(1),y0*theta(1),theta(3))*alpha)/phi; beta_pr=(normpdf(y(1),y0*theta(2),theta(4))*beta)/phi; LogL= LogL+log(phi); %====================================================== for i=2:T % Predicción % Pr(S(t+1)=0jphi(t)) alpha = p1*alpha_pr+(1-p1)*beta_pr; % Pr(S(t+1)=1jphi(t)) beta = (1-q1)*alpha_pr+q1*beta_pr; %====================================================== % Actualización phi=normpdf(y(i),y(i-1)*theta(1),theta(3))*alpha + normpdf(y(i),y(i-1)*theta(2),theta(4))*beta; alpha_pr=(normpdf(y(i),y(i-1)*theta(1),theta(3))*alpha)/phi; beta_pr=(normpdf(y(i),y(i-1)*theta(2),theta(4))*beta)/phi; LogL=LogL+log(phi); end LogL=-LogL;
172CAPíTULO 12 CADENAS DE MARKOV Y MODELOS DE CAMBIOS DE RÉG
end
12.6.3. Filter.m function [Filt_pr,Filt_act]=Filter(theta) % Filtrado de probabilidades no observadas (predicción y actualización). % Se asume una matriz de transición invariante en el tiempo. % Los parámetros theta deben ser previamente estimados (ML). % alpha: Pr(St=0jphi(t-1)) % beta: Pr(St=1jphi(t-1)) % alpha_pr: Pr(St=0jphi(t)) % beta_pr: Pr(St=1jphi(t)) % Fernando Pérez Forero % 02/10/05 global y; y0=randn(1); p1=1/(1+exp(theta(5))); q1=1/(1+exp(theta(6))); P=[p1,(1-p1);(1-q1),q1]; Pi=markov_ss(P); T=length(y); Filt_pr=zeros(T,2); Filt_act=zeros(T,2); %t=1 alpha = p1*Pi(1,1)+(1-p1)*Pi(2,1); beta = (1-q1)*Pi(1,1)+q1*Pi(2,1); Filt_pr(1,1)=alpha; Filt_pr(1,2)=beta; phi=normpdf(y(1),y0*theta(1),theta(3))*alpha+normpdf(y(1),y0*theta(2),theta(4))*beta; alpha_pr=(normpdf(y(1),y0*theta(1),theta(3))*alpha)/phi; beta_pr=(normpdf(y(1),y0*theta(2),theta(4))*beta)/phi; Filt_act(1,1)=alpha_pr; Filt_act(1,2)=beta_pr; %======================================================
12.6 APÉNDICE: PROGRAMAS DE MATLAB
173
for i=2:T % Predicción % Pr(S(t+1)=0jphi(t)) alpha = p1*alpha_pr+(1-p1)*beta_pr; % Pr(S(t+1)=1jphi(t)) beta = (1-q1)*alpha_pr+q1*beta_pr; Filt_pr(i,1)=alpha; Filt_pr(i,2)=beta; %====================================================== % Actualización phi=normpdf(y(i),y(i-1)*theta(1),theta(3))*alpha+ normpdf(y(i),y(i-1)*theta(2),theta(4))*beta; alpha_pr=(normpdf(y(i),y(i-1)*theta(1),theta(3))*alpha)/phi; beta_pr=(normpdf(y(i),y(i-1)*theta(2),theta(4))*beta)/phi; Filt_act(i,1)=alpha_pr; Filt_act(i,2)=beta_pr; end end
12.6.4. Smooth.m function Smooth_m=Smooth(Filt_pr,Filt_act,theta) % Corrección de probabilidades con información de toda la muestra. % La iteración es hacia atrás. % Las probabilidades condicionales al pasado deben ser estimadas % anteriormente (Filter). % alpha_pr: Pr(St=0jphi(T)) % beta_pr: Pr(St=1jphi(T)) % Fernando Pérez Forero % 02/10/05 global y; y0=randn(1); p1=1/(1+exp(theta(5))); q1=1/(1+exp(theta(6))); P=[p1,(1-p1);(1-q1),q1]; T=length(y); Smooth_m=zeros(T,2);
174CAPíTULO 12 CADENAS DE MARKOV Y MODELOS DE CAMBIOS DE RÉG %t=T alpha = Filt_act(T,1); beta = Filt_act(T,2); Smooth_m(T,1)=alpha; Smooth_m(T,2)=beta; for i=T:-1:2 alpha = Smooth_m(i,1); beta = Smooth_m(i,2); alpha_pr=(alpha*p1*Filt_act(i-1,1))/Filt_pr(i,1)+(beta*(1-p1)*Filt_act(i-1,1))/Filt_pr(i,2); beta_pr=(alpha*(1-q1)*Filt_act(i-1,2))/Filt_pr(i,1)+(beta*q1*Filt_act(i-1,2))/Filt_pr(i,2); Smooth_m(i-1,1)=alpha_pr; Smooth_m(i-1,2)=beta_pr; end end
12.6.5. markov_ss.m function [Pi,V,Val]=markov_ss(P) % Cálculo de la distribución estacionaria de probabilidades del estado. % El argumento es la matriz de transición P. % Pi : Distribución estacionaria. % V : Vectores propios de P. % Val: Vector columna con los valores propios de P. % Fernando Pérez Forero % 02/10/05 [V,D]=eig(P); m=length(P); T=eye(m); A=T-P’; clear T; y=ones(1,m); A=[A;y]; clear y; Tr=((A’*A)^(-1))*A’; Pi=Tr(:,m+1); Val=diag(D);
12.6 APÉNDICE: PROGRAMAS DE MATLAB
175
12.6.6. simulate_markov.m1 function [chain,state] = simulate_markov(x,P,pi0,T); % % x = the quantity corresponding to each state, typical element x(i) % % P = Markov transition matrix, typical element p(i,j) i,j=1,...n % % pi0 = probability distribution over initial state %% T = number of periods to simulate %% % % chain = sequence of realizations from the simulation %% Modi…cation of progam by L&S. n = length(x); % % what is the size of the state vector? E = rand(T,1); % % T-vector of draws from independent uniform [0,1] cumsumP = P*triu(ones(size(P))); %% creates a matrix whose rows are the cumulative sums of %% the rows of P % % % % % SET INITIAL STATE USING pi0 E0 = rand(1,1); ppi0 = [0,cumsum(pi0)]; s0 = ((E0<=ppi0(2:n+1)).*(E0>ppi0(1:n)))’; s = s0; %%%%% ITERATE ON THE CHAIN for t=1:T, state(:,t) = s; ppi = [0,s’*cumsumP]; s = ((E(t)<=ppi(2:n+1)).*(E(t)>ppi(1:n)))’; end chain = x’*state;
12.6.7. ejemplo_markov.m % Ejemplo de simulación de una Cadena de Markov % Fernando Pérez Forero % 02/10/05 P=[0.7,0.3;0.2,0.8]; % Matriz estocástica de Transición. s=[0;1];% Valores que toma el estado. 1
Tomado del web site de Chris Edmond: http://pages.stern.nyu.edu/~cedmond/406/
176CAPíTULO 12 CADENAS DE MARKOV Y MODELOS DE CAMBIOS DE RÉG T=50; % Número de períodos a simular. pi0=markov_ss(P); pi0=pi0’;% Distribución inicial(estado estacionario). % Simulación [cadena,estado] = simulate_markov(s,P,pi0,T); …gure(1) plot(cadena,’bo-’) title(’Simulación de una Cadena de Markov’)
12.6.8. tauchen.m2 function [prob,eps,z]=tauchen(N,mu,ro,sig); % Discretizes an AR(1) process into a Markov chain. Determines the optimal grid % and transition matrix. Based on Tauchen (1991). % % y(t)= mu(1-ro) + ro*y(t-1) + u(t) with V(u) = sig^2 % % syntax: % % [prob,eps,z]=tauchen(N,mu,ro,sig) % % N is the number of points on the grid. % prob is the transition matrix % eps contains the cut-o¤ points from - infty to + infty % z are the grid points, i.e. the conditional mean within [eps(i),eps(i+1)]. global mu_ ro_ sigEps_ sig_ eps_ jindx_ if N==1;prob=1; eps=mu;z=mu; else; if ro==0; sigEps=sig; eps=sigEps*norminv((0:N)/N)+mu; eps(1)=-20*sigEps+mu; eps(N+1)=20*sigEps+mu; aux=(eps-mu)/sigEps; aux1=aux(1:end-1); aux2=aux(2:end); 2
Tomado de www.eco.utexas.edu/~cooper/dynprog/tauchen.m
12.6 APÉNDICE: PROGRAMAS DE MATLAB z=N*sigEps*(normpdf(aux1)-normpdf(aux2))+mu; prob=ones(N,N)/N; else; sigEps=sig/sqrt(1-ro^2); eps=sigEps*norminv((0:N)/N)+mu; eps(1)=-20*sigEps+mu; eps(N+1)=20*sigEps+mu; aux=(eps-mu)/sigEps; aux1=aux(1:end-1); aux2=aux(2:end); z=N*sigEps*(normpdf(aux1)-normpdf(aux2))+mu; mu_=mu;ro_=ro;sigEps_=sigEps;eps_=eps;sig_=sig; prob=zeros(N,N); for i=1:N for jindx_=1:N prob(i,jindx_)=quadl(@integ3,eps_(i),eps_(i+1),1e-6)*N; end end end z=z’; eps=eps’; end function F=integ3(u); global mu_ ro_ sigEps_ eps_ jindx_ sig_ aux1=(eps_(jindx_)-mu_*(1-ro_)-ro_*u)/sig_; aux2=(eps_(jindx_+1)-mu_*(1-ro_)-ro_*u)/sig_; F=(normcdf(aux2)-normcdf(aux1)); F=F.*exp(-0.5*(u-mu_).*(u-mu_)/sigEps_^2); pi=4*atan(1); F=F/sqrt(2*pi*sigEps_^2);
12.6.9. consumo.m % Modelado de preferencias vía Cadenas de Markov % Basado en L&S (2000), Ej 1.3 % Fernando Pérez Forero % 02/10/05 c = [1;5];
177
