Predicción del IBEX mediante índice temporal
Francisco Parra
Doctor en Economía
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.
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)
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”
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")
}
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=α+β_1t+β_2t2+u, y se construye el regresor x=α+β_1t+β_2t2 que se representa junto con su periodograma los crecimientos 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 β y α 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")
Bibliografía
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, (https://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.
Predicción del IBEX mediante indice temporal (R-Pub)
Comparación con pronósticos automáticos
Econometria aplicada 