Francisco Parra
Introducción
Utilizando la metodología de la regresión dependiente del tiempo o si se prefiere de la frecuencia (Parra, 2012), se realiza una predicción del crecimiento del IBEX español para el primer trimestre del 2013 al cuarto del 2013, utilizando la serie de crecimientos anuales del IBEX, correspondiente al periodo primer trimestre de 1996 a IV trimestre de 2012. Como predictor se utiliza una tendencia cuadrática que es ajustada utilizando mínimos cuadrados ordinarios.
El análisis dedica unos primeros apartados a exponer la metodología estadística utilizada, y un aparatado en donde se realiza la predicción y lo realmente ocurrido.
Analisis del espectro
Nerlove (1964) y Granger (1969) fueron los primeros investigadores en aplicar el analisis espectral a las series de tiempo en economía. El uso del analisis espectral requiere un cambio en el modo de ver las series económicas, al pasaser de la perspectiva del tiempo al dominio de la frecuencia. El analisis espectral parte de la supsición de que cuanquier serie {Xt}, puede ser transformada en ciclos formados con senos u cosenos:
\[ \begin{equation} X_t=\eta+\sum_{j=1}^N[a_j\cos(2\pi\frac{ft}n)+b_j\sin(2\pi\frac{ft}n)] \end{equation} \] (1)
donde \( \eta \) es la media de la serie, \( a_j \) y \( b_j \) son su amplitud,\( f \) son las frecuencias que del conjunto de las \( n \) observaciones,\( t \) es un indice de tiempo que va de 1 a N, siendo N el numero de periodos para los cuales tenemos observaciones en el conjunto de datos, el cociente \( \frac{ft}n) \) convierte cada valor de \( t \) en escala de tiempo en proporciones de \( 2n \) y rango \( j \) desde \( 1 \) hasta \( n \) siendo \( n=\frac{N}2 \) (es decir, 0,5 ciclos por intervalo de tiempo). Las dinámica de las altas frecuencias (los valores más altos de f) corresponden a los ciclos cortos en tanto que la dinámica de la bajas frecuencias (pequeños valores de f) van a corresponder con los ciclos largos. Si nosotros hacemos que \( \frac{ft}n=w \) la ecuación (1) quedaría, asi :
\[ \begin{equation} X_t=\eta+\sum_{j=1}^N[a_j\cos(\omega_j)+b_j\sin(\omega_j)] \end{equation} \](2)
El analisis espectral puede utilizarse para identificar y cuantificar en procesos aparentemente aperiodicos, sucesiones de cicos de periodo de corto y largo plazo. Una serie dada \( {X_t} \) puede contener diversos ciclos de diferentes frecuencias y amplitudes, y esa combinación de frecuencias y amplitudes de carcter cíclico la hacer aparecer como un serie no periodica e irregular. De hecho la ecuación (2), muestra que cada observación \( t \) de una serie de tiempo, es el resultado sumar los valores en \( t \) que resultan de N ciclos de diferente longitud y amplitud, a los que habría que añadir si cabe un termino de error.
Una manera practica de pasar desde el dominio del tiempo al dominio de la frecuencia es pre-multiplicando los datos originales por una matriz ortogonal, W, sugerida por Harvey (1978), con el elemento (j,t)th :
\[ \begin{equation} w_{jt} = \left\lbrace \begin{array}{ll} \left(\frac{1}T\right) ^\frac{1}2 & \forall j=1\\ \left(\frac{2}T\right) ^\frac{1}2 \cos\left[\frac{\pi j(t-1)}T\right] & \forall j=2,4,6,..\frac{(T-2)}{(T-1)}\\ \left(\frac{2}T\right) ^\frac{1}2 \sin\left[\frac{\pi (j-1)(t-1)}T\right] & \forall j=3,5,7,..\frac{(T-2)}T\\ \left(\frac{1}T\right) ^\frac{1}2 (-1)^{t+1} & \forall j=T \end{array} \right. \end{equation} \](3)
La matriz \( W \) tiene la ventaja de ser ortogonal por lo que \( WW^T=I \).
Función gdf(a) y Función gdt (a)
La función gdf(a) transforma los datos del dominio del tiempo al dominio de la frecuencia pre-multiplicandolos por la matriz ortogonal,\( W \), sugerida por Harvey (1978) y la función gdt (a) transforma los datos del dominio de la frecuencia al dominio del tiempo.
gdf <- function(a) {
a <- matrix(a, nrow = 1)
n <- length(a)
uno <- as.numeric(1:n)
A <- matrix(rep(sqrt(1/n), n), nrow = 1)
for (i in 3:n - 1) {
if (i%%2 == 0) {
A1 <- matrix(sqrt(2/n) * cos(pi * i * (uno - 1)/n), nrow = 1)
A <- rbind(A, A1)
} else {
A2 <- matrix(sqrt(2/n) * sin(pi * (i - 1) * (uno - 1)/n), nrow = 1)
A <- rbind(A, A2)
}
}
AN <- matrix(sqrt(1/n) * (-1)^(uno + 1), nrow = 1)
A <- rbind(A, AN)
A %*% t(a)
}
gdt <- function(a) {
a <- matrix(a, nrow = 1)
n <- length(a)
uno <- as.numeric(1:n)
A <- matrix(rep(sqrt(1/n), n), nrow = 1)
for (i in 3:n - 1) {
if (i%%2 == 0) {
A1 <- matrix(sqrt(2/n) * cos(pi * i * (uno - 1)/n), nrow = 1)
A <- rbind(A, A1)
} else {
A2 <- matrix(sqrt(2/n) * sin(pi * (i - 1) * (uno - 1)/n), nrow = 1)
A <- rbind(A, A2)
}
}
AN <- matrix(sqrt(1/n) * (-1)^(uno + 1), nrow = 1)
A <- rbind(A, AN)
t(A) %*% t(a)
}
Function cdt (a)
Otiene la matriz auxiliar para operaciones con vectores en dominio de tiempo y dominio de la frecuencia, pre-multiplica un vector por la matriz ortogonal, \( W \) y por su transpuesta, Parra F. (2013)
La multiplicación de dos series armónicas de diferente frecuencia:
\[ \begin{equation} [a_j\cos (\omega_j)+b_j\sin (\omega_j)]x [a_i\cos (\omega_i)+b_i\sin (\omega_i)] \end{equation} \]
da como resultado la siguiente suma: \[ \begin{equation} \begin{array}{c} a_ja_i\cos(\omega_j)\cos(\omega_i)+a_jb_i\cos (\omega_j)\sin (\omega_i)\\ +a_ib_j\sin (\omega_j)\cos (\omega_i)b_i\sin (\omega_i)+b_jb_i\sin(\omega_j)\sin(\omega_i) \end{array} \end{equation} \]
considerando las identidades del producto de senos y cosenos, quedaría:
\[ \begin{equation} \begin{array}{c} \frac{a_ja_i+b_jb_i}{2} \cos(\omega_j- \omega_i)+\frac{b_ja_i-b_ja_j}{2}\sin(\omega_j- \omega_i)\\ +\frac{a_ja_i-b_jb_i}{2}\cos(\omega_j+ \omega_i)+\frac{b_ja_i+b_ja_i}{2}\sin(\omega_j+ \omega_i) \end{array} \end{equation} \]
La circularidad de \( \omega \) determina que la serie producto de dos series en \( t \), resulte una nueva serie cuyos coeficientes de Fourier sean una combinación lineal de los coeficientes de Fourier de las series multiplos.
Partiendo de las dos series siguientes:
\[ \begin{equation} \begin{array} {cc} y_t=\eta^y+a_0^y\cos(\omega_0)+b_0^y\sin(\omega_0)+a_1^y\cos(\omega_1)+b_1^y\sin(\omega_1)+ a_2^y\cos(\omega_2)+b_2^y\sin(\omega_2)+a_3^y\cos(\omega_3)\\ x_t=\eta^x+a_0^x\cos(\omega_0)+b_0^x\sin(\omega_0)+a_1^x\cos(\omega_1)+b_1^x\sin(\omega_1)+ a_2^x\cos(\omega_2)+b_2^x\sin(\omega_2)+a_3^x\cos(\omega_3) \end{array} \end{equation} \]
Dada una matriz \( \Theta^{\dot x\dot x} \) de tamaño 8x8 :
\[ \Theta^{\dot x\dot x} = \eta^x I_8+\frac{1}2\left( \begin{array}{cccccccc} 0& a_0^x& b_0^x & a_1^x & b_1^x & a_2^x & b_2^x& 2a_3^x \\ 2a_0^x& a_1^x& b_1^x & a_0^x+a_2^x & b_0^x+b_2^x & a_1^x+2a_3^x & b_1^x& 2a_2^x \\ 2b_0^x& b_1^x&- a_1^x & -b_0^x+b_2^x & a_0^x-a_2^x &- b_1^x &a_1^x- a_3^x &- 2b_2^x \\ 2a_1^x& a_0^x+a_2^x&- b_0^x+b_2^x & 2a_3^x &0 & a_0^x+a_2^x & b_0^x-b_2^x& 2a_1^x \\ 2b_1^x& a_0^x+b_2^x&- b_0^x-a_2^x &0& -2a_3^x & -b_0^x+b_2^x & a_0^x-a_2^x& -2b_1^x \\ 2a_2^x& a_1^x+2a_3^x&- b_1^x & a_0^x+a_2^x &-b_0^x-b_2^x & a_1^x &- b_1^x& 2a_0^x \\ 2b_2^x& b_1^x& a_1^x-2a_3^x & b_0^x-b_2^x &a_0^x-a_2^x & -b_1^x &- a_1^x& -2b_0^x \\ 2a_3^x& a_2^x& -b_2^x & a_1^x &- b_1^x & a_0^x & -b_0^x& 0 \end{array} \right) \] Se demuestra que:
\[ \dot z=\Theta^{\dot x\dot x}\dot y \]
donde \( \dot y = Wy \),\( \dot x = Wx \), y \( \dot z = Wz \).
En el dominio del tiempo:
\[ z_t= x_t y_t=W^T\dot x W^T\dot y=W^T Wx_t W^T\dot y=x_tI_nW^T\dot y \]
\[ W^T\dot z=x_tI_nW^T\dot y \]
\[ \dot z=Wx_tI_nW^T\dot y \]
Entonces:
\[ \Theta^{\dot x\dot x}=W^Tx_tI_nW \]
La matriz cuadrada \( \Theta^{\dot x\dot x} \) puede ser utilizada para obtener los resultados en el dominio de la frecuencia de diversas funciones de series de tiempo . Por ejemplo, si se desea obtener el desarrollo de los coeficientes en fourier de \( z_t=x_t^2 \), entonces:
\( \dot z= Wx_tI_nW^T\dot x \)
En consecuencia, si \( z_t=x_t^n \)
\( \dot z= Wx_t^{n-1}I_nW^T\dot x \)
Si ahora queremos obtener el desarrollo en coeficientes de fourier de \( z_t=\frac{x_t}{y_t} \), entonces:
\( \dot z= W[\frac{1}y_t]I_nW^T\dot x \)
cdf <- function(a) {
a <- matrix(a, nrow = 1)
n <- length(a)
uno <- as.numeric(1:n)
A <- matrix(rep(sqrt(1/n), n), nrow = 1)
for (i in 3:n - 1) {
if (i%%2 == 0) {
A1 <- matrix(sqrt(2/n) * cos(pi * i * (uno - 1)/n), nrow = 1)
A <- rbind(A, A1)
} else {
A2 <- matrix(sqrt(2/n) * sin(pi * (i - 1) * (uno - 1)/n), nrow = 1)
A <- rbind(A, A2)
}
}
AN <- matrix(sqrt(1/n) * (-1)^(uno + 1), nrow = 1)
A <- rbind(A, AN)
I <- diag(c(a))
B <- A %*% I
B %*% t(A)
}
Función periodograma (a)
Calcula y presenta el espectro de la serie “a”
Sea \( a \) un vector n x 1 el modelo transformado en el dominio de la frecuencia esta dado por:
\( \hat a= Wa \)
Denominando \( p_j \) el ordinal del periodograma de \( \hat a \) en la frecuencia \( \lambda_j=2\pi j/n \), y \( \hat a_j \) el j-th elemento de \( \hat a \), entonces
\[ \left\lbrace \begin{array}{ll} p_j=\hat a_{2j}^{2}+\hat a_{2j+1}^{2} & \forall j = 1,...\frac{n-1}{2}\\ p_j=\hat a_{2j}^{2}& \forall j = \frac{n}{2}-1 \end{array} \right . \]
\[ p_0=\hat a_{1}^{2} \]
Entonces el cuadrado del \( \hat a \) puede ser utilizado como un estimador consistente del periodograma de \( a \).
periodograma <- function(a) {
cf <- gdf(a)
n <- length(a)
if (n%%2 == 0) {
m1 <- c(0)
m2 <- c()
for (i in 1:n) {
if (i%%2 == 0)
m1 <- c(m1, cf[i]) else m2 <- c(m2, cf[i])
}
m2 <- c(m2, 0)
frecuencia <- seq(0:(n/2))
frecuencia <- frecuencia - 1
omega <- pi * frecuencia/(n/2)
periodos <- n/frecuencia
densidad <- (m1^2 + m2^2)/(4 * pi)
tabla <- data.frame(omega, frecuencia, periodos, densidad)
tabla$densidad[(n/2 + 1)] <- 2 * tabla$densidad[(n/2 + 1)]
data.frame(tabla[2:(n/2 + 1), ])
} else {
m1 <- c(0)
m2 <- c()
for (i in 1:(n - 1)) {
if (i%%2 == 0)
m1 <- c(m1, cf[i]) else m2 <- c(m2, cf[i])
}
m2 <- c(m2, cf[n])
frecuencia <- seq(0:((n - 1)/2))
frecuencia <- frecuencia - 1
omega <- pi * frecuencia/(n/2)
periodos <- n/frecuencia
densidad <- (m1^2 + m2^2)/(4 * pi)
tabla <- data.frame(omega, frecuencia, periodos, densidad)
data.frame(tabla[2:((n + 1)/2), ])
}
}
Función gperiodogrma (a)
Presenta gráficamente el espectro de la variabe a
gperiodograma <- function(a) {
tabla <- periodograma(a)
plot(tabla$frecuencia, tabla$densidad, main = "Espectro", ylab = "densidad",
xlab = "frecuencia", type = "l", col = "#ff0000")
}
Regresión band spectrum
Realiza la regresión band spectrum del vector de datos a con el vector de datos b para las frecuencias que figuran en el vector c.
Hannan (1963) fue quien propuso la regresión en dominio de la frecuencia (regresión band spectrum). Engle (1974), demostró que dicha regresión no alteraba los supuestos básicos de la regresión clásica, cuyos estimadores eran Estimadores Lineales Insesgados y Optimos (ELIO).
\[ \begin{equation} y=X\beta+u \end{equation} \](4)
donde \( X \) es una matriz \( n x k \) de observaciones de \( k \) variables independientes, \( \beta \) es un vecto \( k x I \) de parámetros, \( y \) es un vecto \( n x 1 \) de observaciones de la variable dependiente, y \( u \) en un vector \( n x I \) de pertubacciones de media cero y varianza constante, \( \sigma^2 \).
El modelo puese expresarse en el dominio de la frecuencia aplicando una transformación lineal a las variables dependiente e independientes,por ejemplo, premultiplicando todas las variables por la matriz ortogonal \( W \). La técnica de la “regresión band spectrum”,consiste en realizar el analisis de regresión en el dominio de la frecuencia pero omitinedo determinadas oscilaciones periodicas. Con este procedimiento pueden tratarse problemas derivados de la estacionalidad de las series o de autocorrelación en los residuos. Engle (1974) muesta que si los residuos están correlacionados serialmente y son generados por un procieso estacionario estocastico, la regresión en el dominio de la frecuencia es el estimador asintóticamente más eficiente de \( \beta \).
La transforamción de la ecuación (4) del dominio del tiempo al dominio de la frecuencia en Engle (1974), se bas en la matriz W, cuyo elemento \( (t, s) \) esta dado por:
\( w_{ts}=\frac{1}{\sqrt n} e^{i\lambda_t s} \),\( s= 0,1,...,n-1 \)
donde \( \lambda_t = 2\pi \frac{t}n \), \( t=0,1,...,n-1 \), y \( i=\sqrt{-1} \).
Premultiplicando las observaciones de (4) por W, obtenemos: \[ \begin{equation} \dot y=\dot X\beta+\dot u \end{equation} \](5)
donde \( \dot y = Wy \),\( \dot X = WX \), y \( \dot u = Wu \).
Si el vector de las perturbaciones en (4) cumple las hipoteis clásicas del modelo de regresión: \( E[u] = 0 \) y \( E[uu']=\sigma^2 I_n \). entonces el vector de perturbaciones transformado al dominio de la frecuencia, \( \dot u \), tendrá las mims propiedades. Por otro lado, dado que la matriz W es ortogonal., \( WW^{T}= I \), entonces \( W^T \) sería la transpuesta de la completa conjugada de W. De forma que las observaciones del modelo (5) acaban conteniendo el mismo tipo de información que las observaciones del modelo inicialmente planteado.
Si aplicamos MCO a (5) , dadas las propiedades de \( \dot u \), obtendríamos el mejor estimador lineal insesgando (ELIO) de \( \dot \beta \). El estimador obtenido sería de hecho identico al estimador MCO de (4).
Estimar (5) manteniendo unicamente determinadas frecuencias, puede llevarse a cabo omitiendo las observaciones correspondientes a las restantes frecuencias, si bien, dado que las variables en (5) son complejas, Engle (1974) sugiere la transformada inversa de Fourier para recomponer el modelo estimado en términos de tiempo.
Definiendo una matriz de tamaño \( A \) de tamaño n x n de ceros excepto en las posiciones de la diagonal principal correspondientes a las frecuencias que queremos incluir en la regresión y premultiplicando \( \dot y \) y \( \dot X \) por \( A \) eleminamos determindas observaciones y las reemplazamos por ceros para realizar la regresión band spectrum. Devolver al dominio del tiempo estas observaciones requiere:
\[ \begin{equation} y^* = W^{T}A\dot y = W^{T}AWy \\ x^* = W^{T}A\dot x = W^{T}AWx \end{equation} \](6)
Al regresar \( y^* \) sobre \( x^* \) obtenemos un \( \beta \) identico al estimador que obtendríamos al estimar por MCO \( \dot y \) frente a \( \dot x \).
Cuando se realiza la regresión band spectrum de esta mnera, ocurre un problema asociado a los grados de libertad de la regresión de \( y^* \) sobre \( x^* \) que asumen los programas estadisticos convencionales, \( n - k \), en vez de los grados de libertad reales que serían unicamente \( n'- k \), donde \( n' \) es el numero de frecuencias incluidas en la regresión band spectrum.
Si la regresión espectral va a ser usada para obtener un estimador asintóticamente eficiente de \( \beta \) en presencia de autocorrelación en el termino de error, la matriz \( A \) ha de ser reemplazada por otra matriz diagonal, \( V \). En dicha diagonal principal ha de incluirse el estimador de \( \int_u^{1/2}(\lambda) \), donde \( \int_u (\lambda) \) es el la transformación del termino de error obtenido en (4) al dominio de la frecuencia \( \lambda \). Puede utilizarse un programa convencional para obtener \( \beta \) haciendo una transformación análoga a (6); ver Engle and Gardner [1976]. Sin embargo, si el procedimiento va a ser iterativo, lo que podría llevar a una mejora en las propiededes de las muestra pequeñas, la transformada inversa de fourier debería emplearse en cada iteración.
La regresión con parámetros dependiente del tiempoo de la frecuencia .
Consideramos el modelo de regresión siguiente:
\[ \begin{equation} y_t=\beta_tx_t+u_t \end{equation} \](7)
donde \( x_t \) es un vector T x 1 de observaciones de las variable independiente, \( \beta_t \) es un vector de T x 1 del parametro \( \beta \) , \( y_t \) es un vector de T x 1 observaciones de la variable depenendiente, y \( u_t \) es un vector de errores distribuidos con media cero y varianza constante.
Asumiendo que las series, \( y_t \),\( x_t \),\( \beta_t \) and \( u_t \), pueden ser transformadas en el dominio de la frecuencia:
\[ y_t=\eta^y+\sum_{j=1}^N[a^y_j\cos(\omega_j)+b^y_j\sin(\omega_j) \]
\[ x_t=\eta^x+\sum_{j=1}^N[a^y_j\cos(\omega_j)+b^y_j\sin(\omega_j)] \]
\[ \beta_t=\eta^\beta+\sum_{j=1}^N[a^\beta_j\cos(\omega_j)+b^\beta_j\sin(\omega_j)] \]
Otenemos dichas series pre-multiplicando (7) por \( W \)
\[ \dot y=\dot x\dot\beta+\dot u \]
donde \( \dot y = Wy \),\( \dot x = Wx \), \( \dot \beta = W\beta \) y \( \dot u = Wu \)
El sistema (7) puede reescribirse como (ver anexo):
\[ \begin{equation} \dot y=Wx_tI_nW^T\dot \beta + WI_nW^T\dot u \end{equation} \] (8)
Si denominamos \( \dot e=WI_nW^T\dot u \), podrían buscarse los \( \dot \beta \) que minimizaran la suma cuadrática de los errores \( e_t=W^T\dot e \).
\( min \sum_{t=1}^T [W^T\dot e]^2 \) (9)
Utilizando los mínimos cuadrados ordinarios, encontramos una solución a (9), construyendo una matriz de regresores \( X \) cuya primera columna sería el vector de tamaño \( T \) \( (1,0,0,...)_T \), la segunda columna sería la primera fila de la matriz \( \Theta^{\dot x\dot x} \) y las columnas siguientes, corresponderían las filas de \( \Theta^{\dot x\dot x} \) correspondientes a las frecuencias de senos o cosenos que queremos regresar.
Los coeficietes de la solución MCO: \( \dot \beta = (x^Tx)^{-1}x^T\dot y \) serían: \( \dot \beta_0 \) el asociado a la constante, \( \dot \beta_1 \) el asociado a la pendiente, \( \dot \beta_2 ... \) los asociados a las frecuencias de senos y cosenos elegidas.
La solución mínimocuadrática del modelo de regresión simple \( y_t=\hat \beta_0+\hat \beta_1 x_t+\hat u_t \), sería \( W^T \dot\beta = W^T[(x^Tx)^{-1}x^T\dot y] \).
donde \( x \) es una matriz \( T \) x \( 2 \) cuya primera columna es el vector \( (1,0,0,...)_{T} \), la segunda columna sería la primera fila de la matriz \( \Theta^{\dot x\dot x} \), nótese que dicha columna viene a ser el vector transpuesto de coeficientes de fourier de \( x_t \) dividido entre dos y encabezado por el valor medio de \( x_t \), de manera que el producto matricial una vez desarrollado viene a ser la covarianza de \( x_t \) e \( y_t \) y la varianza de \( x_t \).
Si se trata de una regresión multiple \( y_t=\hat \beta_o+\hat \beta_1 x_t+ \hat \beta_2 z_t + \hat u_t \) la matriz \( x \) incluiría una tercera columna con la primera fila de la matriz \( \Theta^{\dot z\dot z} \).
Estimación del crecimiento del PIB mediante indice
Partimos de los crecimientos interanuales del IBEX español en indices de volumen del periodo I trimestre de 1996 a IV trimestre de 2013, y generamos un indice temporal \( t \) con 67 observaciones. Se realiza la regresión lineal entre \( y = \hat\alpha + \hat\beta_1 t + \hat\beta_2 t^2 + \hat u \), y se construye el regresor \( x = \hat\alpha + \hat\beta_1 t + \hat\beta_2 t^2 \) que se representa junto con su periodograma los crecimeintos interanuales del IBEX.
y <- read.csv("ibex.csv", header = FALSE, sep = ";", dec = ",")
t <- seq(1, 68)
t2 <- t^2
y <- y$V1[1:68]
index <- lm(y ~ t + t2)
x <- index$fitted
plot(ts(y), pch = 19, col = "blue")
points(ts(x), pch = 19, col = "red")
gperiodograma(y)
gperiodograma(x)
El periodograma del IBEX muestra que la mayor parte de su varianza está explicada por las 10 frecuencias más altas, se plantea un modelo de regresión dependiente en el tiempo basado en los 30 primeros coeficientes de la transformación en el dominio de la frecuencia del regresor \( x \), incluyendo constante y pendiente que como vemos como aproxima con bastante precisión los crecimientos interanuales del IBEX.
z <- c(rep(1, 68))
cy <- gdf(y)
cx <- cdf(x)
cz <- cdf(z)
c1 <- matrix(cz[1, ], nrow = 1)
c2 <- matrix(cx[1:30, ], nrow = 30)
X <- rbind(c1, c2)
B <- solve(X %*% t(X)) %*% (X %*% cy)
Y <- t(X) %*% B
fitted <- gdt(Y)
plot(ts((y)), pch = 19, col = "blue")
points(ts((fitted)), pch = 19, col = "red")
Buscando un mejor modelo regresion paso a paso
Se utiliza la función stepAIC del paquete MASS para seleccionar el mejor modelo entre los posibles considerando las frecuencias asociadas a la constante, pendiente y los 30 primeros coeficientes de fourier en la transformación del ibex en el dominio de la frecuencia. Se representa la gráfica temporal del crecimiento del IBEX y la estimación en el dominio del tiempo de la ecuación seleccionada, el periodograma acumulado de los residuos del modelo y el diagrama de dispersión entre el crecimiento del IBEX ,\( Y \), y el indice temporal,\( X \), en donde se incluye tambien el resultado de la regresión seleccionada.
rownames(X) <- c(paste("X", 1:31, sep = "_"))
X <- data.frame(t(X))
library(MASS)
fit <- lm(cy ~ 0 + X_1 + X_2 + X_3 + X_4 + X_5 + X_6 + X_7 + X_8 + X_9 + X_10 +
X_11 + X_12 + X_13 + X_15 + X_16 + X_17 + X_18 + X_19 + X_20 + X_21 + X_22 +
X_23 + X_25 + X_26 + X_27 + X_28 + X_29 + X_30 + X_31, X)
step <- stepAIC(fit, direction = "both")
## Start: AIC=297
## cy ~ 0 + X_1 + X_2 + X_3 + X_4 + X_5 + X_6 + X_7 + X_8 + X_9 +
## X_10 + X_11 + X_12 + X_13 + X_15 + X_16 + X_17 + X_18 + X_19 +
## X_20 + X_21 + X_22 + X_23 + X_25 + X_26 + X_27 + X_28 + X_29 +
## X_30 + X_31
##
## Df Sum of Sq RSS AIC
## - X_19 1 20 2303 296
## - X_22 1 47 2331 296
## - X_29 1 52 2336 296
## - X_15 1 55 2338 297
## - X_23 1 66 2350 297
## <none> 2284 297
## - X_7 1 101 2385 298
## - X_31 1 118 2402 298
## - X_3 1 132 2416 299
## - X_20 1 194 2478 300
## - X_28 1 205 2488 301
## - X_26 1 224 2508 301
## - X_30 1 240 2524 302
## - X_6 1 312 2596 304
## - X_18 1 329 2613 304
## - X_16 1 383 2667 305
## - X_10 1 427 2711 307
## - X_17 1 467 2751 308
## - X_25 1 470 2754 308
## - X_27 1 481 2764 308
## - X_13 1 520 2804 309
## - X_12 1 596 2880 311
## - X_21 1 834 3118 316
## - X_11 1 943 3227 318
## - X_8 1 1129 3413 322
## - X_9 1 1159 3443 323
## - X_5 1 1478 3762 329
## - X_1 1 2418 4702 344
## - X_2 1 2769 5053 349
## - X_4 1 3232 5516 355
##
## Step: AIC=295.5
## cy ~ X_1 + X_2 + X_3 + X_4 + X_5 + X_6 + X_7 + X_8 + X_9 + X_10 +
## X_11 + X_12 + X_13 + X_15 + X_16 + X_17 + X_18 + X_20 + X_21 +
## X_22 + X_23 + X_25 + X_26 + X_27 + X_28 + X_29 + X_30 + X_31 -
## 1
##
## Df Sum of Sq RSS AIC
## - X_15 1 35 2339 295
## - X_29 1 41 2345 295
## <none> 2303 296
## - X_3 1 113 2417 297
## - X_31 1 114 2417 297
## + X_19 1 20 2284 297
## - X_23 1 148 2451 298
## - X_20 1 174 2478 298
## - X_28 1 186 2489 299
## - X_7 1 225 2528 300
## - X_30 1 235 2538 300
## - X_22 1 262 2565 301
## - X_26 1 326 2630 303
## - X_25 1 451 2754 306
## - X_6 1 508 2811 307
## - X_27 1 523 2826 307
## - X_16 1 606 2910 309
## - X_17 1 653 2956 311
## - X_18 1 672 2976 311
## - X_13 1 748 3051 313
## - X_12 1 769 3072 313
## - X_10 1 773 3076 313
## - X_21 1 827 3130 314
## - X_8 1 1130 3433 321
## - X_9 1 1294 3597 324
## - X_11 1 1538 3841 328
## - X_5 1 1818 4121 333
## - X_1 1 2423 4726 342
## - X_2 1 3862 6165 360
## - X_4 1 4256 6559 365
##
## Step: AIC=294.6
## cy ~ X_1 + X_2 + X_3 + X_4 + X_5 + X_6 + X_7 + X_8 + X_9 + X_10 +
## X_11 + X_12 + X_13 + X_16 + X_17 + X_18 + X_20 + X_21 + X_22 +
## X_23 + X_25 + X_26 + X_27 + X_28 + X_29 + X_30 + X_31 - 1
##
## Df Sum of Sq RSS AIC
## - X_29 1 51 2389 294
## <none> 2339 295
## - X_3 1 97 2435 295
## + X_15 1 35 2303 296
## - X_31 1 116 2455 296
## - X_23 1 121 2459 296
## + X_19 1 0 2338 297
## - X_28 1 213 2552 298
## - X_22 1 227 2565 299
## - X_30 1 241 2579 299
## - X_20 1 269 2607 300
## - X_26 1 301 2640 301
## - X_7 1 375 2713 303
## - X_25 1 498 2836 306
## - X_27 1 513 2852 306
## - X_16 1 595 2934 308
## - X_6 1 652 2991 309
## - X_12 1 781 3119 312
## - X_17 1 782 3121 312
## - X_13 1 880 3219 314
## - X_21 1 964 3303 316
## - X_18 1 1135 3474 319
## - X_8 1 1178 3516 320
## - X_10 1 1244 3583 322
## - X_9 1 1266 3605 322
## - X_5 1 2189 4528 337
## - X_1 1 2493 4832 342
## - X_11 1 2523 4861 342
## - X_2 1 4665 7004 367
## - X_4 1 5199 7537 372
##
## Step: AIC=294
## cy ~ X_1 + X_2 + X_3 + X_4 + X_5 + X_6 + X_7 + X_8 + X_9 + X_10 +
## X_11 + X_12 + X_13 + X_16 + X_17 + X_18 + X_20 + X_21 + X_22 +
## X_23 + X_25 + X_26 + X_27 + X_28 + X_30 + X_31 - 1
##
## Df Sum of Sq RSS AIC
## - X_31 1 69 2458 294
## <none> 2389 294
## + X_29 1 51 2339 295
## + X_15 1 45 2345 295
## - X_3 1 100 2489 295
## + X_19 1 4 2385 296
## - X_30 1 193 2582 297
## - X_28 1 199 2589 297
## - X_23 1 221 2611 298
## - X_20 1 233 2622 298
## - X_22 1 317 2706 301
## - X_7 1 385 2774 302
## - X_25 1 473 2862 304
## - X_26 1 518 2907 305
## - X_16 1 573 2963 307
## - X_27 1 628 3017 308
## - X_6 1 653 3042 308
## - X_17 1 748 3137 311
## - X_12 1 772 3161 311
## - X_13 1 859 3248 313
## - X_21 1 916 3305 314
## - X_18 1 1191 3580 320
## - X_8 1 1195 3584 320
## - X_9 1 1270 3659 321
## - X_10 1 1272 3661 321
## - X_5 1 2226 4615 337
## - X_1 1 2534 4924 341
## - X_11 1 2561 4950 342
## - X_2 1 4690 7079 366
## - X_4 1 5260 7650 371
##
## Step: AIC=294
## cy ~ X_1 + X_2 + X_3 + X_4 + X_5 + X_6 + X_7 + X_8 + X_9 + X_10 +
## X_11 + X_12 + X_13 + X_16 + X_17 + X_18 + X_20 + X_21 + X_22 +
## X_23 + X_25 + X_26 + X_27 + X_28 + X_30 - 1
##
## Df Sum of Sq RSS AIC
## <none> 2458 294
## + X_31 1 69 2389 294
## - X_3 1 89 2547 294
## + X_15 1 40 2418 295
## - X_30 1 125 2583 295
## + X_29 1 4 2455 296
## + X_19 1 2 2457 296
## - X_28 1 199 2658 297
## - X_23 1 215 2674 298
## - X_20 1 232 2691 298
## - X_22 1 307 2765 300
## - X_7 1 380 2839 302
## - X_25 1 469 2927 304
## - X_26 1 475 2933 304
## - X_27 1 569 3027 306
## - X_16 1 580 3039 306
## - X_6 1 669 3128 308
## - X_17 1 739 3197 310
## - X_12 1 766 3224 310
## - X_13 1 871 3329 313
## - X_21 1 935 3393 314
## - X_8 1 1163 3622 318
## - X_18 1 1177 3636 319
## - X_9 1 1244 3702 320
## - X_10 1 1251 3709 320
## - X_5 1 2185 4643 335
## - X_1 1 2477 4936 339
## - X_11 1 2535 4993 340
## - X_2 1 4685 7144 365
## - X_4 1 5208 7667 369
step$anova # display resultss
## Stepwise Model Path
## Analysis of Deviance Table
##
## Initial Model:
## cy ~ 0 + X_1 + X_2 + X_3 + X_4 + X_5 + X_6 + X_7 + X_8 + X_9 +
## X_10 + X_11 + X_12 + X_13 + X_15 + X_16 + X_17 + X_18 + X_19 +
## X_20 + X_21 + X_22 + X_23 + X_25 + X_26 + X_27 + X_28 + X_29 +
## X_30 + X_31
##
## Final Model:
## cy ~ X_1 + X_2 + X_3 + X_4 + X_5 + X_6 + X_7 + X_8 + X_9 + X_10 +
## X_11 + X_12 + X_13 + X_16 + X_17 + X_18 + X_20 + X_21 + X_22 +
## X_23 + X_25 + X_26 + X_27 + X_28 + X_30 - 1
##
##
## Step Df Deviance Resid. Df Resid. Dev AIC
## 1 39 2284 297.0
## 2 - X_19 1 19.60 40 2303 295.5
## 3 - X_15 1 35.21 41 2339 294.6
## 4 - X_29 1 50.62 42 2389 294.0
## 5 - X_31 1 69.18 43 2458 294.0
step$coefficients
## X_1 X_2 X_3 X_4 X_5 X_6 X_7 X_8 X_9
## 163.911 22.024 3.372 -31.865 -15.808 10.212 6.962 8.625 8.870
## X_10 X_11 X_12 X_13 X_16 X_17 X_18 X_20 X_21
## -10.343 -8.540 -7.142 -6.475 4.496 4.658 -3.128 -2.821 -4.402
## X_22 X_23 X_25 X_26 X_27 X_28 X_30
## -2.463 -2.804 -2.734 3.722 3.446 1.699 -1.195
fitted <- gdt(step$fitted.values)
plot(ts((y)), pch = 19, col = "blue")
points(ts(fitted), pch = 19, col = "red")
cpgram(ts(step$resid))
plot((x), (y), pch = 19, col = "blue")
points((x), (fitted), pch = 19, col = "red")
Los residuos del modelo tienen apariencia aleatoria, según la función cpgram del paquete MASS.
Los parametros beta obtenidos una vez se recuperan en el dominio del tiempo, constiuyen en realidad un ajuste del cociente entre el crecimiento del IBEX y del índice construido, tal y como puede verse en la siguiente gráfica:
cero <- data.frame(c(rep(0, 69)))
rownames(cero) <- paste("X", 1:69, sep = "_")
coef <- data.frame(beta = step$coefficients)
coef <- data.frame(coef, id = rownames(coef))
cero <- data.frame(cero, id = rownames(cero))
beta1 <- merge(coef, cero, by = "id", all.y = TRUE)
beta2 <- data.frame(id = beta1$id, beta = beta1$beta, id2 = as.numeric(substr(beta1$id,
3, 4)))
beta2[is.na(beta2)] <- 0
Orden_B = order(beta2$id2)
BaseOrdenada = beta2[Orden_B, ]
beta <- c(BaseOrdenada$beta[2:69])
alpha <- c(BaseOrdenada$beta[1:1], rep(0, 67))
beta_t <- gdt(beta)
alpha_t <- gdt(alpha)
fitted1 <- alpha_t + beta_t * (x)
plot(((y) - alpha_t)/(x), pch = 19, col = "blue")
lines(beta_t, pch = 19, col = "red")
Se representa el periodograma de los coeficientes temporales beta obtenidos, en la selección del modelo paso a paso.
gperiodograma(beta_t)
Prediccion
Al ser los coeficientes beta circulares, la predicción se realiza con los coeficientes correspondientes a las primeras observaciones, con ellos se obtienen los nuevos coeficientes de fourier para los parametros \( \dot\beta \) y \( \dot\alpha \) considerando ahora un \( N=72 \):
beta2_t <- c(beta_t, beta_t[1, 1], beta_t[2, 1], beta_t[3, 1], beta_t[4, 1])
beta2_f <- gdf(beta2_t)
plot(ts(beta2_f))
points(ts(beta), pch = 19, col = "red")
alpha2_t <- c(alpha_t, alpha_t[1, 1], alpha_t[2, 1], alpha_t[3, 1], alpha_t[4,
1])
alpha2_f <- gdf(alpha2_t)
plot(ts(alpha2_f))
points(ts(alpha), pch = 19, col = "red")
Por ultimo se estiman de nuevo el modelo, se reprentan sus resltados en grafico de serie temporal y en diagrama de dispersión X-Y.
t <- seq(1, 72)
t2 <- t^2
B2 <- c(alpha2_f[1, 1], beta2_f)
x2 <- predict(index, data.frame(t, t2))
z2 <- c(rep(1, 72))
cx2 <- cdf(x2)
cz2 <- cdf(z2)
c1 <- matrix(cz2[1, ], nrow = 1)
X2 <- rbind(c1, cx2)
Y2 <- t(X2) %*% B2
fitted2 <- gdt(Y2)
fitted21 <- t(alpha2_t) + t(beta2_t) * t(x2)
y <- read.csv("ibex.csv", header = FALSE, sep = ";", dec = ",")
plot(ts(y), pch = 19, col = "blue")
points(ts(fitted2), pch = 19, col = "red")
points(ts(fitted), pch = 19, col = "green")
plot(ts(gdt(beta2_f)))
points(ts(beta_t), pch = 19, col = "red")
Bibliography
Engle, Robert F. (1974), Band Spectrum Regression,International Economic Review 15,1-11.
Hannan, E.J. (1963), Regression for Time Series, in Rosenblatt, M. (ed.), Time Series Analysis, New York, John Wiley.
Harvey, A.C. (1978), Linear Regression in the Frequency Domain, International Economic Review, 19, 507-512.
Parra, F. (2013), Regresión con paramétros dependientes del tiempo, (http://econometria.wordpress.com/2013/07/29/estimacion-con-parametros-dependientes-del-tiempo/)
Ripley, B. Venables, B.; Bates, D.M.; Hornik, K. Gebhardt, G. and Firth, D. (2013): “Package MASS”. http://www.stats.ox.ac.uk/pub/MASS4/
Wilson, P.J. and Perry, L.J. (2004). Forecasting Australian Unemployment Rates Using Spectral Analysis Australian Jurnal of Labour Economics, vol 7,no 4, December 2004, pp 459-480.