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 en euros el primer trimestre del 2013 al cuarto del 2013, a partir de la serie correspondiente al periodo primer trimestre de 1996 a IV trimestre de 2012.
La predicción se realia en varias etapas:
a) Primera etapa: Se realiza una regresión MCO del IBEX con un indice temporal \( t \), para obtener una tendencia, primero Se introduce un tendencia logaritmica \( y_t=\hat\alpha + \hat\beta_1t + \hat\beta_2t^2 + \hat u_t \).
b) Segunda etapa: Se elige el porcentaje de la varianza del cociente entre la serie y la tendencia que se quiere explicar, y se plantea una regresión band-spectrum utilizando la transformación de \( y_t \) en el dominio de la frecuencia (función gdf) y la circularización de la tendencia cuadrática (función cdt (a)), el resultado de esta regresión son uno parametros dependientes del tiempo \( \hat \beta_t \) que al ser multiplicados por la tendencia del crecimiento IBEX la aproximan a la evolución real del crecimiento del IBEX.
c) Tercera etapa: mediante una regresión paso a paso, se selecciona el mejor modelo eliminando los regresores \( \dot \beta \) irrelevantes, eligiendose un modelo restringido que es el que se utiliza para predecir.
4) Cuarta etapa: Predicción. Al ser los coeficientes \( \hat \beta_t \) circulares, la predicción para \( t+k \) se realiza a través del producto de los \( k \) primeros coeficientes y los \( t+k \) valores de la tendencia ajustada.
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. Los pronosticos se comparan con los pronosticos de los modelos automáticos ARIMA y ETS del paquete “forecast”.
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)
Obtiene 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 tiempo/ de la frecuencia .
Consideramos ahora 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 IBEX mediante tendencia polinómica
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 68 observaciones. Se realiza la regresión lineal entre \( y_t=\hat\alpha + \hat\beta_1t + \hat\beta_2t^2 + \hat u_t \), y se construye el regresor. Se representa el IBEX junto a su tendencia, tambien el gráfico del periodograma del cociente \( \frac{y_t}{\hat \alpha + \hat \beta_1t + \hat \beta_2t} \).
y <- read.csv("ibex.csv", header = FALSE, sep = ";", dec = ",")
n <- 68
t <- seq(1, n)
t2 <- t^2
y <- y$V1[1:n]
index <- lm(y ~ t + t2)
x <- index$fitted
plot(ts(y), pch = 19, col = "blue")
points(ts(x), pch = 19, col = "red")
u <- y/index$fitted
plot(ts(u), pch = 19, col = "blue")
per <- periodograma(u)
gperiodograma(u)
p <- as.numeric(per$densidad)
n2 <- length(p)
s <- p[1]
t <- 1:n2
for (i in 2:n2) {
s1 <- p[i] + s[(i - 1)]
s <- c(s, s1)
s2 <- s/s[n2]
}
S <- data.frame(frecuencia = per$frecuencia, s2)
sel <- subset(S, s2 > 0.9)
m <- c(sel[1, 1] * 2)
Se plantea un modelo de regresión en bandas de frecuencia basado utilizando los coeficientes de la transformación en el dominio de la frecuencia del regresor \( x \),vinculados al 90% de la varianza del cociente entre la serie y la tendencia, incluyendo constante y pendiente que como vemos aproxima con bastante precisión los crecimientos interanuales del IBEX.
z <- c(rep(1, n))
cy <- gdf(y)
cx <- cdf(x)
cz <- cdf(z)
c1 <- matrix(cz[1, ], nrow = 1)
c2 <- matrix(cx[1:m, ], nrow = m)
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")
Se utiliza de nuevo 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.Se representa el periodograma acumulado de los residuos del modelo seleccionado con la función cpgram del paquete MASS.
X <- matrix(t(X), nrow = n)
library(MASS)
fit <- lm(cy ~ 0 + X)
step <- stepAIC(fit, direction = "both")
## Start: AIC=132.7
## cy ~ 0 + X
##
## Df Sum of Sq RSS AIC
## <none> 75 133
## - X 63 48221 48296 446
step$anova # display resultss
## Stepwise Model Path
## Analysis of Deviance Table
##
## Initial Model:
## cy ~ 0 + X
##
## Final Model:
## cy ~ 0 + X
##
##
## Step Df Deviance Resid. Df Resid. Dev AIC
## 1 5 75 132.7
step$coefficients
## X1 X2 X3 X4 X5 X6 X7
## 223.31639 21.69958 8.92754 -36.56147 -21.19847 6.36230 7.08559
## X8 X9 X10 X11 X12 X13 X14
## 13.90172 13.59987 -13.63670 -13.01231 -9.94105 -6.25639 4.47286
## X15 X16 X17 X18 X19 X20 X21
## 3.78273 1.75570 1.44361 -4.58782 1.22458 -1.00427 -4.34985
## X22 X23 X24 X25 X26 X27 X28
## -4.73622 -3.71098 1.23236 -1.77892 3.09851 3.33436 1.84449
## X29 X30 X31 X32 X33 X34 X35
## 1.45864 -2.80523 -2.41038 -1.39994 -0.45026 1.43976 0.17159
## X36 X37 X38 X39 X40 X41 X42
## -1.12597 -1.71101 0.33790 0.09186 1.38382 0.30032 0.03822
## X43 X44 X45 X46 X47 X48 X49
## -0.47201 -0.18888 -0.27906 -0.46602 0.07258 0.10484 1.23424
## X50 X51 X52 X53 X54 X55 X56
## 1.66188 1.96232 -3.11026 -3.10199 0.26399 1.01643 1.48075
## X57 X58 X59 X60 X61 X62 X63
## -0.88773 -1.71308 -2.14211 -0.08084 0.39862 1.12231 -0.32521
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 parámetros \( \dot \beta \) obtenidos en la regresión en bandas de frecuencia, se recuperan en el dominio del tiempo, dado que se estima la relación \( \hat y_t=\hat \alpha +\hat \beta_t x_t \), los parámetros \( \hat \beta_t \) pueden representarse en terminos temporales junto a \( \frac{y_t-\hat \alpha}{x_t} \), lo que viene a ser el cociente entre el crecimiento del IBEX y el índice construido. Se representa dicha relación.
cero <- data.frame(c(rep(0, (n + 1))))
rownames(cero) <- paste("X", 1:(n + 1), 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,
2, 3)))
beta2[is.na(beta2)] <- 0
Orden_B = order(beta2$id2)
BaseOrdenada = beta2[Orden_B, ]
beta <- c(BaseOrdenada$beta[2:(n + 1)])
alpha <- c(BaseOrdenada$beta[1:1], rep(0, (n - 1)))
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")
Estimaciones automáticas con el paquete forecast
library(forecast)
## Warning: package 'forecast' was built under R version 3.0.3
## Loading required package: zoo
## Warning: package 'zoo' was built under R version 3.0.3
##
## Attaching package: 'zoo'
##
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
##
## Loading required package: timeDate
## Warning: package 'timeDate' was built under R version 3.0.3
## This is forecast 5.3
# Automated forecasting using an exponential model
fit1 <- ets(ts(y, start = 1, frequency = 4))
# Automated forecasting using an ARIMA model
fit2 <- auto.arima(ts(y, start = 1, frequency = 4))
pred1 <- forecast(fit1, 4)
pred2 <- forecast(fit2, 4)
plot(ts(y, start = 1, frequency = 4), pch = 19, col = "blue")
points((fit1$fitted), pch = 19, col = "red")
# Automated forecasting using an ARIMA model
plot(ts(y, start = 1, frequency = 4), pch = 19, col = "blue")
points(fit2$x - fit2$residual, pch = 19, col = "red")
Predicción
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 \). Se representan los \( \hat \beta \) en término de dominio y frecuencia, junto a los antes estimados.
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)
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)
data.frame(alpha2_t, alpha2_f, beta2_t, beta2_f)
## alpha2_t alpha2_f beta2_t beta2_f
## 1 27.08 2.298e+02 -0.10940 21.04410
## 2 27.08 -1.599e-14 -0.07825 -0.33739
## 3 27.08 4.552e-15 -0.31939 -36.44658
## 4 27.08 -1.599e-14 0.13284 -17.53644
## 5 27.08 -7.661e-15 0.38708 16.50432
## 6 27.08 -5.329e-15 0.82715 13.13905
## 7 27.08 -8.882e-16 1.48987 5.11439
## 8 27.08 4.263e-14 0.68868 -4.55757
## 9 27.08 2.109e-14 1.56948 -17.97120
## 10 27.08 -3.997e-14 1.17782 -15.13790
## 11 27.08 3.908e-14 0.12430 10.21263
## 12 27.08 -8.882e-16 0.42368 5.96655
## 13 27.08 9.770e-15 -0.63317 6.16538
## 14 27.08 -6.839e-14 -1.15347 2.17845
## 15 27.08 5.373e-14 -0.73854 -5.17756
## 16 27.08 4.174e-14 -0.65713 -4.05220
## 17 27.08 2.842e-14 -0.27381 0.16928
## 18 27.08 -8.571e-14 -0.82026 1.82121
## 19 27.08 -1.243e-14 -0.54545 1.99374
## 20 27.08 2.398e-14 -1.75075 0.66496
## 21 27.08 -4.885e-14 -2.39343 4.63900
## 22 27.08 0.000e+00 -2.29688 5.05480
## 23 27.08 6.084e-14 -3.19795 -2.51528
## 24 27.08 -8.882e-16 -2.66141 0.09834
## 25 27.08 2.620e-14 -2.82318 -3.50870
## 26 27.08 -8.726e-14 -2.99543 -3.39994
## 27 27.08 1.696e-13 -3.83638 -1.33973
## 28 27.08 7.572e-14 -3.81735 0.29297
## 29 27.08 -8.615e-14 -4.76214 3.41971
## 30 27.08 -7.905e-14 -3.73202 2.99787
## 31 27.08 -7.194e-14 -1.29683 -0.95180
## 32 27.08 2.331e-15 -0.84580 -1.23701
## 33 27.08 -8.615e-14 0.92749 -2.32745
## 34 27.08 -4.996e-16 -0.46269 -0.63593
## 35 27.08 -1.066e-14 -2.02901 1.14949
## 36 27.08 -1.148e-13 -1.22151 0.20176
## 37 27.08 5.951e-14 -1.81948 -1.12943
## 38 27.08 9.281e-14 -1.35196 -1.50413
## 39 27.08 -7.638e-14 1.19350 0.72801
## 40 27.08 1.333e-13 -0.86510 0.77884
## 41 27.08 1.368e-13 -0.28673 1.16438
## 42 27.08 -2.442e-15 -1.36575 0.23263
## 43 27.08 4.352e-14 -3.33855 -0.33459
## 44 27.08 2.887e-15 4.36948 -0.46801
## 45 27.08 1.288e-13 3.14979 0.03544
## 46 27.08 5.329e-14 6.54284 -0.18533
## 47 27.08 3.197e-14 14.61108 -0.23332
## 48 27.08 1.652e-13 17.56931 0.59065
## 49 27.08 1.941e-13 24.35998 -0.21289
## 50 27.08 4.174e-14 20.36728 1.11618
## 51 27.08 9.592e-14 19.22773 -0.66726
## 52 27.08 1.332e-14 22.71877 -2.70445
## 53 27.08 1.932e-13 18.88398 -0.51427
## 54 27.08 7.194e-14 14.14543 2.04832
## 55 27.08 3.109e-13 6.86769 3.37535
## 56 27.08 7.994e-15 0.18152 0.64074
## 57 27.08 -6.306e-14 -0.80864 -1.91461
## 58 27.08 -1.137e-13 3.83453 -1.01219
## 59 27.08 -4.663e-14 5.75043 1.62905
## 60 27.08 -2.043e-14 6.01404 1.33262
## 61 27.08 -1.865e-13 4.05531 -0.59916
## 62 27.08 -1.901e-13 2.33018 -0.74169
## 63 27.08 -6.062e-14 5.42056 -2.08784
## 64 27.08 9.770e-14 4.59663 -0.51189
## 65 27.08 1.732e-13 5.81196 0.61347
## 66 27.08 -5.684e-14 6.48357 -0.41434
## 67 27.08 3.553e-15 4.85149 -0.21445
## 68 27.08 -5.596e-14 3.14164 -0.12663
## 69 27.08 -1.231e-13 -0.10940 -0.06620
## 70 27.08 3.446e-13 -0.07825 -0.08942
## 71 27.08 5.843e-13 -0.31939 -0.03024
## 72 27.08 0.000e+00 0.13284 -0.05697
plot(ts(gdt(beta2_f)))
points(ts(beta_t), pch = 19, col = "red")
plot(ts(beta2_f))
points(ts(gdf(beta_t)), pch = 19, col = "red")
Por ultimo se estiman de nuevo el modelo, y se reprentan sus resultados en grafico de serie temporal.
t <- c(seq(1, (n + 4)))
t2 <- t^2
B2 <- c(alpha2_f[1, 1], beta2_f)
x2 <- predict(index, data.frame(t, t2))
z2 <- c(rep(1, (n + 4)))
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(t(fitted21)), pch = 19, col = "red")
points(ts(fitted), pch = 19, col = "green")
lines(ts(x2), pch = 19, col = "black")
Se comparan los pronosticos con los obtenidos con los pronosticados por el ajuste automático ARIMA y ETS del paquete “forecast”.
pred1
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 18 Q1 -7.372 -24.46 9.714 -33.50 18.76
## 18 Q2 -7.372 -31.53 16.790 -44.32 29.58
## 18 Q3 -7.372 -36.96 22.220 -52.63 37.88
## 18 Q4 -7.372 -41.54 26.797 -59.63 44.89
pred2
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 18 Q1 -9.8514 -23.80 4.096 -31.18 11.48
## 18 Q2 3.7903 -15.93 23.515 -26.38 33.96
## 18 Q3 0.1177 -24.04 24.275 -36.83 37.06
## 18 Q4 -13.8058 -41.70 14.089 -56.47 28.85
fitted21[68:72]
## [1] -5.226 28.248 27.944 30.721 25.520
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.
Hyndman, R. ; Khandakar Y. (2008): Automatic Time Series Forecasting: The forecast Package for R. Journas of statistical Software. July,2008,Volumen 37,Issue 3.(http://cran.r-project.org/web/packages/forecast/index.html)
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.