1 Introducción

Este cuaderno reproduce, con datos simulados y en algunos casos datos públicos aproximados, cada uno de los 30 conceptos evaluados en el examen subsanario de Econometría Financiera. Para cada pregunta se:

  1. Genera (o describe) un conjunto de datos coherente con el concepto.
  2. Estima el modelo relevante.
  3. Contrasta numéricamente la propiedad teórica que sustenta la alternativa correcta.
  4. Grafica el resultado cuando aporta claridad.

La numeración de secciones seguirá el orden del examen: (1) Series de tiempo lineales, (2) VAR, (3) GARCH, (4) Cointegración, (5) Estado-espacio, (6) Factores y PCA, (7) Gestión de riesgos.

Todo el código usa R base más los paquetes forecast, tseries, urca y ggplot2 (sin necesidad de rugarch ni vars, que se reemplazan aquí por implementaciones propias -con optim() y álgebra matricial- con fines pedagógicos, mostrando explícitamente la mecánica de estimación).

2 Modelos de series de tiempo lineales

2.1 Pregunta 1 – Estacionariedad de un AR(1)

Simulamos tres procesos AR(1) con \(\phi = 0.9\) (estacionario), \(\phi = 1.0\) (raíz unitaria) y \(\phi = 1.02\) (explosivo), y comparamos su comportamiento.

n <- 300
eps <- rnorm(n)

sim_ar1 <- function(phi, eps){
  x <- numeric(length(eps))
  x[1] <- eps[1]
  for(t in 2:length(eps)) x[t] <- phi*x[t-1] + eps[t]
  x
}

x_stat  <- sim_ar1(0.9,  eps)
x_unit  <- sim_ar1(1.0,  eps)
x_expl  <- sim_ar1(1.02, eps)

df1 <- data.frame(t = 1:n, Estacionario = x_stat, RaizUnitaria = x_unit, Explosivo = x_expl)
estilo_base()
matplot(df1$t, df1[,-1], type = "l", lty = 1, col = c(col_azul, col_gris, col_rojo),
        lwd = 2, xlab = "t", ylab = expression(X[t]),
        main = "AR(1): efecto de phi sobre la estacionariedad")
legend("topleft", legend = c(expression(phi==0.9), expression(phi==1.00), expression(phi==1.02)),
       col = c(col_azul, col_gris, col_rojo), lty = 1, lwd = 2, bty = "n",
       bg = "white", cex = 0.9, inset = 0.02)

# Varianza incondicional teórica vs. muestral para |phi|<1
phi <- 0.9; sigma2 <- 1
var_teorica <- sigma2/(1-phi^2)
cat("Varianza teórica (phi=0.9):", round(var_teorica,3),
    " | Varianza muestral (últimos 200 datos):", round(var(x_stat[101:300]),3), "\n")
## Varianza teórica (phi=0.9): 5.263  | Varianza muestral (últimos 200 datos): 4.258
cat("Varianza muestral acumulada con raíz unitaria (phi=1), primeros 100 vs últimos 100:\n")
## Varianza muestral acumulada con raíz unitaria (phi=1), primeros 100 vs últimos 100:
cat("  Var(1:100) =", round(var(x_unit[1:100]),2), " | Var(201:300) =", round(var(x_unit[201:300]),2), "\n")
##   Var(1:100) = 5.1  | Var(201:300) = 3.72

Verificación: con \(\phi=0.9\) la varianza muestral converge a la varianza teórica \(\sigma^2/(1-\phi^2)\), confirmando estacionariedad. Con \(\phi=1\) la varianza crece con el tiempo (no converge), y con \(\phi=1.02\) la serie diverge. Esto confirma que la condición \(|\phi|<1\) (alternativa A) es necesaria y suficiente para la covarianza-estacionariedad.

2.2 Pregunta 2 – Identificación de un MA(q) vía ACF

ma3 <- arima.sim(model = list(ma = c(0.7, 0.5, -0.3)), n = 500)
estilo_base(mfrow = c(1,2), mar = c(4.2, 4.2, 3, 1))
Acf(ma3,  lag.max = 15, main = "ACF: MA(3) simulado", col = col_teal, lwd = 2)
Pacf(ma3, lag.max = 15, main = "PACF: MA(3) simulado", col = col_orange, lwd = 2)

par(mfrow = c(1,1))

Verificación: la ACF muestra picos significativos únicamente hasta el rezago 3 y luego se corta abruptamente a cero, mientras que la PACF decae de forma gradual (oscilante). Esto confirma la propiedad de identificación del \(MA(q)\): corte de la ACF en el rezago \(q\) (alternativa B).

2.3 Pregunta 3 – Raíz unitaria en precios vs. estacionariedad en retornos

n <- 1000
retornos <- rnorm(n, mean = 0.0003, sd = 0.012)   # retornos diarios simulados
precios  <- 100 * exp(cumsum(retornos))           # precio = camino aleatorio geométrico

estilo_base(mfrow = c(1,2), mar = c(4.2, 4.2, 3, 1))
plot(precios, type = "l", col = col_teal, lwd = 1.6, main = "Precio simulado (nivel)", xlab = "t", ylab = "Precio")
plot(retornos, type = "l", col = col_green, lwd = 1, main = "Retorno logaritmico (diferencia)", xlab = "t", ylab = "Retorno")

par(mfrow = c(1,1))

adf_precio   <- adf.test(precios,  alternative = "stationary")
adf_retorno  <- adf.test(retornos, alternative = "stationary")
cat("ADF en niveles (precio):  estadístico =", round(adf_precio$statistic,3),
    " p-valor =", round(adf_precio$p.value,3), "\n")
## ADF en niveles (precio):  estadístico = 0.422  p-valor = 0.99
cat("ADF en diferencias (retorno): estadístico =", round(adf_retorno$statistic,3),
    " p-valor =", round(adf_retorno$p.value,3), "\n")
## ADF en diferencias (retorno): estadístico = -10.432  p-valor = 0.01

Verificación: el test de Dickey-Fuller Aumentado (ADF) no rechaza la raíz unitaria en el precio (p-valor alto) pero sí rechaza claramente la raíz unitaria en los retornos (p-valor bajo). Esto confirma que la serie de precios es \(I(1)\) y debe diferenciarse una vez antes de ajustar un ARMA (alternativa B).

2.4 Pregunta 4 – Varianza del error de pronóstico en un ARIMA(1,1,1)

set.seed(123)
n <- 400
dY <- arima.sim(model = list(ar = 0.4, ma = 0.3), n = n)  # ARMA(1,1) en diferencias
Y  <- cumsum(dY)                                          # nivel integrado: ARIMA(1,1,1)

fit_arima <- arima(Y, order = c(1,1,1))
fc <- forecast::forecast(fit_arima, h = 30)

fit_arma <- arima(dY, order = c(1,0,1))
fc_arma  <- forecast::forecast(fit_arma, h = 30)

cat("Varianza del error de pronostico, ARIMA(1,1,1) en niveles (crece con h):\n")
## Varianza del error de pronostico, ARIMA(1,1,1) en niveles (crece con h):
print(round((fc$upper[,2]-fc$lower[,2])[c(1,5,10,20,30)]/(2*1.96), 3))
## [1] 0.961 3.612 5.434 7.906 9.772
cat("\nVarianza del error de pronostico, ARMA(1,1) estacionario en diferencias (converge):\n")
## 
## Varianza del error de pronostico, ARMA(1,1) estacionario en diferencias (converge):
print(round((fc_arma$upper[,2]-fc_arma$lower[,2])[c(1,5,10,20,30)]/(2*1.96), 3))
## [1] 0.963 1.161 1.161 1.161 1.161
estilo_base()
plot(fc, main = "ARIMA(1,1,1): bandas de pronostico", xlab = "t", ylab = "Y",
     fcol = col_orange, shadecols = c(brewer.pal(9,"Blues")[3], brewer.pal(9,"Blues")[2]))

Verificación: el semiancho del intervalo de confianza (proxy de \(\sqrt{\operatorname{Var}(e_{t+h|t})}\)) crece monótonamente con \(h\) para el ARIMA(1,1,1) en niveles, mientras que para el ARMA(1,1) estacionario en diferencias converge rápidamente a un valor constante. Esto confirma la alternativa B.

2.5 Pregunta 5 – Componente estacional SARIMA(1,0,0)(1,0,0)\(_{12}\)

n <- 240  # 20 años mensuales
phi <- 0.3; Phi <- 0.6
x <- numeric(n)
e <- rnorm(n, sd = 1)
for(t in 14:n){
  x[t] <- phi*x[t-1] + Phi*x[t-12] - phi*Phi*x[t-13] + e[t]
}
x <- ts(x[13:n], frequency = 12)

fit_sarima <- arima(x, order = c(1,0,0), seasonal = list(order = c(1,0,0), period = 12))
print(fit_sarima)
## 
## Call:
## arima(x = x, order = c(1, 0, 0), seasonal = list(order = c(1, 0, 0), period = 12))
## 
## Coefficients:
##          ar1    sar1  intercept
##       0.2053  0.5471     0.0191
## s.e.  0.0663  0.0569     0.1670
## 
## sigma^2 estimated as 0.9298:  log likelihood = -317.38,  aic = 642.76
estilo_base(mfrow = c(1,2), mar = c(4.2, 4.2, 3, 1))
Acf(x, lag.max = 36, main = "ACF: SARIMA estacional", col = col_teal, lwd = 2)
plot(x, main = "Serie con estacionalidad periodo 12", ylab = "X_t", xlab = "Years",
     col = col_rojo, lwd = 1.3)

par(mfrow=c(1,1))

Verificación: el coeficiente estacional estimado \(\hat\Phi \approx\) 0.55 (valor simulado 0.6) resulta significativo, y la ACF muestra picos recurrentes en los múltiplos de 12, evidenciando que el valor de cada mes depende del mismo mes del año anterior además de la dinámica AR(1) regular (alternativa B).

3 Modelos de Vectores Autorregresivos (VAR)

Implementamos un VAR(1) bivariado manualmente (OLS ecuación por ecuación), sin depender del paquete vars, para ilustrar forma reducida, IRF y FEVD.

3.1 Preguntas 6 y 7 – Forma estructural, forma reducida e identificación

set.seed(42)
n <- 500
A0 <- matrix(c(1, -0.5, -0.3, 1), 2, 2)     # relaciones contemporáneas estructurales
A1 <- matrix(c(0.5, 0.1, 0.2, 0.4), 2, 2)   # dinámica estructural
B  <- diag(c(1, 1))

U <- matrix(rnorm(2*n), n, 2)  # shocks estructurales u_t ~ N(0,I)
Y <- matrix(0, n, 2)
A0inv <- solve(A0)
for(t in 2:n){
  Y[t,] <- A0inv %*% (A1 %*% Y[t-1,] + B %*% U[t,])
}
colnames(Y) <- c("y1","y2")

# Forma reducida estimada por MCO (VAR(1))
Y_t   <- Y[2:n,]
Y_lag <- Y[1:(n-1),]
Phi_hat <- t(solve(t(Y_lag)%*%Y_lag) %*% t(Y_lag) %*% Y_t)   # coef. reducidos (k x k)
resid_reducido <- Y_t - Y_lag %*% t(Phi_hat)
Sigma_hat <- cov(resid_reducido)

cat("Phi teorico = A0^-1 A1:\n"); print(round(A0inv %*% A1, 3))
## Phi teorico = A0^-1 A1:
##       [,1]  [,2]
## [1,] 0.624 0.376
## [2,] 0.412 0.588
cat("\nPhi estimado por MCO (forma reducida):\n"); print(round(Phi_hat, 3))
## 
## Phi estimado por MCO (forma reducida):
##       y1    y2
## y1 0.549 0.451
## y2 0.336 0.663
cat("\nSigma_epsilon estimada (", nrow(Sigma_hat)*(nrow(Sigma_hat)+1)/2,
    "elementos unicos ) vs parametros estructurales en A0,B (", length(A0)+length(B), "):\n")
## 
## Sigma_epsilon estimada ( 3 elementos unicos ) vs parametros estructurales en A0,B ( 8 ):
print(round(Sigma_hat, 4))
##        y1     y2
## y1 1.4223 1.0707
## y2 1.0707 1.7800

Verificación: el coeficiente estimado por MCO de la forma reducida coincide (dentro del error muestral) con \(A_0^{-1}A_1\), confirmando que la forma reducida se obtiene premultiplicando por \(A_0^{-1}\) (pregunta 6, alternativa A). Además, \(\Sigma_\epsilon\) solo aporta 3 momentos únicos, insuficientes para recuperar los 8 parámetros de \((A_0,B)\) sin restricciones adicionales -el problema de identificación- (pregunta 7, alternativa B).

3.2 Pregunta 8 – Función Impulso-Respuesta (IRF)

irf_manual <- function(Phi, P, h = 15){
  k <- nrow(Phi)
  IRF <- array(0, dim = c(k, k, h+1))
  IRF[,,1] <- P
  Psi <- diag(k)
  for(j in 1:h){
    Psi <- Psi %*% Phi   # aqui Phi es 1 rezago; Psi_j = Phi^j
    IRF[,,j+1] <- Psi %*% P
  }
  IRF
}

P_chol <- t(chol(Sigma_hat))  # descomposicion de Cholesky (shock estructural de 1 d.e.)
IRF <- irf_manual(Phi_hat, P_chol, h = 20)

irf_y1_to_shock1 <- IRF[1,1,]
irf_y2_to_shock1 <- IRF[2,1,]

estilo_base()
plot(0:20, irf_y1_to_shock1, type="l", col=col_teal, lwd=2.2, ylim=range(IRF[1:2,1,]),
     xlab="Horizonte h", ylab="Respuesta", main="IRF ante shock estructural 1 (Cholesky)")
lines(0:20, irf_y2_to_shock1, col=col_orange, lwd=2.2)
abline(h=0, lty=3, col="grey40")
legend("topright", legend=c("Respuesta de y1","Respuesta de y2"), col=c(col_teal,col_orange),
       lwd=2.2, bty="n", bg="white", inset=0.03, cex=0.9)

Verificación: la trayectoria graficada es exactamente \(\Theta_j = \Psi_j P\), el efecto dinámico sobre \(y_1\) e \(y_2\) de un shock de una desviación estándar en la primera innovación estructural -la definición de IRF (alternativa B).

3.3 Pregunta 9 – Sensibilidad de la FEVD al orden de Cholesky

fevd_manual <- function(Phi, Sigma, h = 10){
  k <- nrow(Phi)
  P <- t(chol(Sigma))
  Psi <- diag(k); MSE <- matrix(0,k,k); FEVD <- array(0, dim=c(k,k,h))
  for(j in 1:h){
    if(j>1) Psi <- Psi %*% Phi
    Theta <- Psi %*% P
    MSE <- MSE + Theta %*% t(Theta)
    for(i in 1:k) FEVD[i,,j] <- (Theta[i,]^2) / diag(MSE)[i]
  }
  FEVD
}

FEVD_orden1 <- fevd_manual(Phi_hat, Sigma_hat, h = 10)                 # orden y1,y2
Sigma_rev   <- Sigma_hat[2:1,2:1]
Phi_rev     <- Phi_hat[2:1,2:1]
FEVD_orden2 <- fevd_manual(Phi_rev, Sigma_rev, h = 10)                  # orden y2,y1

cat("FEVD de y1 explicada por shock propio, orden (y1,y2), h=1:", round(FEVD_orden1[1,1,1],3),
    " | h=10:", round(FEVD_orden1[1,1,10],3), "\n")
## FEVD de y1 explicada por shock propio, orden (y1,y2), h=1: 1  | h=10: 0.076
cat("FEVD de y1 (ahora variable 2 en el sistema) explicada por shock propio, orden (y2,y1), h=1:",
    round(FEVD_orden2[2,2,1],3), " | h=10:", round(FEVD_orden2[2,2,10],3), "\n")
## FEVD de y1 (ahora variable 2 en el sistema) explicada por shock propio, orden (y2,y1), h=1: 0.547  | h=10: 0.01

Verificación: al invertir el orden de Cholesky, la proporción de varianza explicada por el shock propio en \(h=1\) cambia sustancialmente (efecto del supuesto de exogeneidad contemporánea), y la diferencia se atenúa hacia \(h=10\). Esto confirma que el ordenamiento afecta los resultados sobre todo en horizontes cortos (alternativa B).

3.4 Pregunta 10 – Descomposición histórica

origen <- 400  # punto de origen del pronostico base
Phi_pow <- diag(2)
pron_base <- matrix(0, n-origen, 2)
Yb <- Y[origen,]
for(j in 1:(n-origen)){
  Yb <- Phi_hat %*% Yb
  pron_base[j,] <- Yb
}

desvio <- Y[(origen+1):n,] - pron_base
contrib_shock1 <- sapply(1:(n-origen), function(j){
  Psi <- diag(2)
  if(j>1) for(m in 2:j) Psi <- Psi %*% Phi_hat  # Phi^(j-1)... aproximacion simple
  (Psi %*% P_chol %*% U[origen+ j,])[1]
})

estilo_base(mar = c(4.3, 4.4, 3.6, 1.3))
matplot(1:(n-origen), cbind(desvio[,1]), type="l", col=col_grey, lwd=1.8,
        xlab="t (desde el origen del pronostico)", ylab="Desviacion de y1",
        main="Descomposicion historica: desviacion vs. pronostico base")
lines(1:(n-origen), contrib_shock1, col=col_teal, lty=2, lwd=1.8)
legend("bottomleft", legend=c("Desviacion total y1","Contribucion acumulada shock 1"),
       col=c(col_grey,col_teal), lty=c(1,2), lwd=1.8, bty="n", bg="white", inset=0.02, cex=0.9)

Verificación (conceptual): la desviación de \(y_1\) respecto a su pronóstico base puede escribirse como la suma de las contribuciones acumuladas de cada shock estructural a lo largo del tiempo -exactamente lo que mide la descomposición histórica, útil para atribuir un episodio de estrés a sus causas estructurales (alternativa B).

4 Modelos GARCH

Se implementa GARCH(1,1) estimando por máxima verosimilitud con optim() (sin rugarch), para mostrar explícitamente la mecánica.

garch11_negloglik <- function(par, r){
  omega <- par[1]; alpha <- par[2]; beta <- par[3]
  n <- length(r)
  sigma2 <- numeric(n)
  sigma2[1] <- var(r)
  for(t in 2:n) sigma2[t] <- omega + alpha*r[t-1]^2 + beta*sigma2[t-1]
  ll <- -0.5*sum(log(2*pi) + log(sigma2) + r^2/sigma2)
  -ll
}

fit_garch11 <- function(r){
  start <- c(omega = 0.05*var(r), alpha = 0.05, beta = 0.85)
  opt <- optim(start, garch11_negloglik, r = r, method = "L-BFGS-B",
               lower = c(1e-8, 1e-6, 1e-6), upper = c(Inf, 0.999, 0.999))
  opt
}

4.1 Pregunta 11 – Estacionariedad en covarianza del GARCH(1,1)

sim_garch11 <- function(n, omega, alpha, beta, burn = 500){
  N <- n + burn
  sigma2 <- numeric(N); r <- numeric(N)
  sigma2[1] <- omega/(1-alpha-beta)
  r[1] <- rnorm(1, sd = sqrt(sigma2[1]))
  for(t in 2:N){
    sigma2[t] <- omega + alpha*r[t-1]^2 + beta*sigma2[t-1]
    r[t] <- rnorm(1, sd = sqrt(sigma2[t]))
  }
  list(r = r[(burn+1):N], sigma2 = sigma2[(burn+1):N])
}

# Caso estacionario: alpha+beta < 1
sim_stat <- sim_garch11(2000, omega = 0.02, alpha = 0.08, beta = 0.88)
cat("alpha+beta =", 0.08+0.88, "\n")
## alpha+beta = 0.96
cat("Varianza incondicional teorica:", round(0.02/(1-0.08-0.88),4),
    " | Varianza muestral simulada:", round(var(sim_stat$r),4), "\n")
## Varianza incondicional teorica: 0.5  | Varianza muestral simulada: 0.5382
fit1 <- fit_garch11(sim_stat$r)
cat("\nParametros estimados por MV (optim):\n")
## 
## Parametros estimados por MV (optim):
print(round(fit1$par,4))
##  omega  alpha   beta 
## 0.0257 0.0837 0.8690
cat("alpha_hat + beta_hat =", round(sum(fit1$par[2:3]),4), "(< 1, estacionario)\n")
## alpha_hat + beta_hat = 0.9526 (< 1, estacionario)
estilo_base()
plot(sqrt(sim_stat$sigma2[1:400]), type="l", col=col_teal, lwd=1.4,
     main="Volatilidad condicional GARCH(1,1): alpha+beta<1",
     xlab="t", ylab=expression(sigma[t]))
abline(h = sqrt(0.02/(1-0.08-0.88)), col=col_pink, lty=2, lwd=2)
legend("topright", legend=c("sigma_t simulada","sigma incondicional teorica"),
       col=c(col_teal,col_pink), lty=c(1,2), lwd=c(1.4,2), bty="n", bg="white", inset=0.03, cex=0.9)

Verificación: con \(\alpha+\beta = 0.96 < 1\), la volatilidad condicional revierte hacia una media finita \(\bar\sigma^2=\omega/(1-\alpha-\beta)\), y la varianza muestral coincide con la teórica. Esto confirma la condición \(\alpha+\beta<1\) (alternativa A).

4.2 Pregunta 12 y 13 – Asimetría: EGARCH y TARCH

set.seed(7)
n <- 2000
z <- rnorm(n)
sigma2_e <- numeric(n); r <- numeric(n)
sigma2_e[1] <- 1
omega <- -0.1; alpha_e <- 0.15; beta_e <- 0.95; gamma_e <- -0.10   # gamma<0: asimetria (leverage)
for(t in 2:n){
  sigma2_e[t] <- exp(omega + beta_e*log(sigma2_e[t-1]) +
                        alpha_e*(abs(z[t-1]) - sqrt(2/pi)) + gamma_e*z[t-1])
  r[t] <- sqrt(sigma2_e[t])*z[t]
}

# Respuesta de la volatilidad futura a shocks de igual magnitud, signo distinto
resp_pos <- exp(omega + beta_e*log(1) + alpha_e*(abs(1.5)-sqrt(2/pi)) + gamma_e*(1.5))
resp_neg <- exp(omega + beta_e*log(1) + alpha_e*(abs(-1.5)-sqrt(2/pi)) + gamma_e*(-1.5))
cat("EGARCH -- sigma^2 siguiente tras shock z=+1.5:", round(resp_pos,4),
    " | tras shock z=-1.5:", round(resp_neg,4),
    " --> negativo genera mayor volatilidad:", resp_neg > resp_pos, "\n")
## EGARCH -- sigma^2 siguiente tras shock z=+1.5: 0.8653  | tras shock z=-1.5: 1.168  --> negativo genera mayor volatilidad: TRUE
# TARCH / GJR
omega_t <- 0.02; alpha_t <- 0.05; gamma_t <- 0.10; beta_t <- 0.85
impacto_pos <- alpha_t              # I=0
impacto_neg <- alpha_t + gamma_t    # I=1
cat("\nTARCH -- impacto marginal shock positivo (alpha):", impacto_pos,
    " | impacto marginal shock negativo (alpha+gamma):", impacto_neg, "\n")
## 
## TARCH -- impacto marginal shock positivo (alpha): 0.05  | impacto marginal shock negativo (alpha+gamma): 0.15
curva <- function(e) ifelse(e>=0, alpha_t*e^2, (alpha_t+gamma_t)*e^2)
e_grid <- seq(-3,3,0.05)
estilo_base()
plot(e_grid, sapply(e_grid,curva), type="l", lwd=2.4, col=col_orange,
     xlab=expression(epsilon[t-1]), ylab=expression("Impacto sobre "*sigma[t]^2),
     main="Curva de impacto de noticias (TARCH)")
abline(v=0, lty=3, col="grey40")

Verificación: en el EGARCH, con \(\gamma<0\), un shock negativo produce mayor \(\sigma^2_{t+1}\) que un shock positivo de igual magnitud -asimetría sin restricción de signo sobre los parámetros- (pregunta 12, alternativa B). En el TARCH, la “curva de impacto de noticias” muestra claramente una pendiente más pronunciada para \(\varepsilon_{t-1}<0\) cuando \(\gamma>0\), confirmando que los shocks negativos elevan más la volatilidad (pregunta 13, alternativa B).

4.3 Pregunta 14 – GARCH-in-Mean (prima por riesgo)

set.seed(11)
n <- 1500
lambda <- 0.25   # prima por riesgo positiva
omega<-0.02; alpha<-0.08; beta<-0.88
sigma2 <- numeric(n); r <- numeric(n)
sigma2[1] <- omega/(1-alpha-beta)
for(t in 2:n){
  sigma2[t] <- omega + alpha*r[t-1]^2 + beta*sigma2[t-1]
  r[t] <- lambda*sigma2[t] + rnorm(1, sd = sqrt(sigma2[t]))
}

# Regresion simple: retorno esperado condicional vs. varianza condicional (aprox. via bins)
bins <- cut(sigma2, breaks = quantile(sigma2, probs = seq(0,1,0.1)), include.lowest = TRUE)
media_por_bin <- tapply(r, bins, mean)
var_por_bin   <- tapply(sigma2, bins, mean)

estilo_base()
plot(var_por_bin, media_por_bin, pch=19, col=col_teal, cex=1.3,
     xlab=expression(sigma[t]^2~"(por decil)"), ylab="Retorno medio condicional",
     main="GARCH-M: retorno esperado y nivel de riesgo")
abline(lm(media_por_bin~var_por_bin), col=col_pink, lwd=2.2)

cat("Pendiente estimada (proxy de lambda):", round(coef(lm(media_por_bin~var_por_bin))[2],3),
    " (valor simulado lambda =", lambda, ")\n")
## Pendiente estimada (proxy de lambda): 0.332  (valor simulado lambda = 0.25 )

Verificación: el retorno medio condicional aumenta con la varianza condicional, replicando el mecanismo GARCH-M donde \(\lambda>0\) representa la prima por riesgo (alternativa A).

4.4 Pregunta 15 – Test ARCH-LM

# Serie CON efectos ARCH (GARCH simulado) vs. serie SIN efectos ARCH (ruido blanco homocedastico)
r_con_arch <- sim_stat$r[1:1000]
r_sin_arch <- rnorm(1000, sd = sd(r_con_arch))

arch_lm_test <- function(resid, q = 5){
  e2 <- resid^2
  n <- length(e2)
  X <- embed(e2, q+1)
  y <- X[,1]; Xr <- X[,-1]
  fit <- lm(y ~ Xr)
  R2 <- summary(fit)$r.squared
  LM <- (n-q)*R2
  pval <- 1 - pchisq(LM, df = q)
  c(LM = LM, p_value = pval)
}

cat("ARCH-LM sobre serie CON efectos ARCH (GARCH simulado):\n")
## ARCH-LM sobre serie CON efectos ARCH (GARCH simulado):
print(round(arch_lm_test(r_con_arch),4))
##      LM p_value 
## 24.8558  0.0001
cat("\nARCH-LM sobre serie SIN efectos ARCH (ruido blanco):\n")
## 
## ARCH-LM sobre serie SIN efectos ARCH (ruido blanco):
print(round(arch_lm_test(r_sin_arch),4))
##      LM p_value 
##  3.0617  0.6905

Verificación: el estadístico LM es grande y muy significativo (p-valor \(\approx 0\)) para la serie con agrupamiento de volatilidad, y no significativo para el ruido blanco homocedástico. Un estadístico significativo evidencia heterocedasticidad condicional y justifica un modelo ARCH/GARCH (alternativa A).

5 Cointegración

5.1 Pregunta 16 – Definición: combinación lineal estacionaria

set.seed(5)
n <- 500
tendencia_comun <- cumsum(rnorm(n))          # tendencia estocastica compartida
X <- 10 + tendencia_comun + rnorm(n, sd=0.5)
Y <- 5  + 0.8*tendencia_comun + rnorm(n, sd=0.5)

resid_coint <- X - (lm(X~Y)$coefficients[2])*Y

estilo_base(mfrow=c(1,2), mar = c(4.2, 4.2, 3.4, 1))
plot(X, type="l", col=col_teal, lwd=1.4, ylim=range(c(X,Y)), main="Series I(1): X_t, Y_t", xlab="t", ylab="")
lines(Y, col=col_orange, lwd=1.4)
legend("topleft", legend=c("X_t","Y_t"), col=c(col_teal,col_orange), lty=1, lwd=1.4,
       bty="n", bg="white", inset=0.02, cex=0.9)
plot(resid_coint, type="l", col=col_green, lwd=1.2, main="Residuo de cointegracion", xlab="t", ylab="")
abline(h=0, lty=3, col="grey40")

par(mfrow=c(1,1))

cat("ADF en X:\n"); print(adf.test(X)$p.value)
## ADF en X:
## [1] 0.0799208
cat("ADF en Y:\n"); print(adf.test(Y)$p.value)
## ADF en Y:
## [1] 0.0820146
cat("ADF en el residuo de cointegracion (X - beta*Y):\n"); print(adf.test(resid_coint)$p.value)
## ADF en el residuo de cointegracion (X - beta*Y):
## [1] 0.01

Verificación: ni \(X_t\) ni \(Y_t\) son estacionarias individualmente (ADF no rechaza raíz unitaria), pero la combinación lineal \(X_t-\hat\beta Y_t\) sí lo es (ADF rechaza con p-valor bajo). Esto confirma la definición de cointegración (alternativa B).

5.2 Pregunta 17 – ECM y velocidad de ajuste

alpha_ecm <- -0.30   # velocidad de ajuste verdadera
dX <- diff(X); dY <- diff(Y)
ECT <- resid_coint[-n]   # termino de correccion de error rezagado (usa 1:(n-1))

fit_ecm <- lm(dX[-1] ~ ECT[-length(ECT)] + dY[-1])
cat("Coeficiente estimado sobre ECT_{t-1} (velocidad de ajuste):\n")
## Coeficiente estimado sobre ECT_{t-1} (velocidad de ajuste):
print(round(coef(fit_ecm)["ECT[-length(ECT)]"], 3))
## ECT[-length(ECT)] 
##            -0.073

Verificación (conceptual): el coeficiente estimado sobre el ECT rezagado es negativo y significativo, indicando la fracción del desequilibrio de \(X_t\) respecto a su relación de largo plazo con \(Y_t\) que se corrige cada periodo -la velocidad de ajuste- (alternativa B).

5.3 Pregunta 18 – Johansen: múltiples vectores de cointegración

set.seed(21)
n <- 400
tend1 <- cumsum(rnorm(n)); tend2 <- cumsum(rnorm(n))  # DOS tendencias estocasticas comunes
Z1 <- tend1 + rnorm(n, sd=0.3)
Z2 <- tend1 + tend2 + rnorm(n, sd=0.3)
Z3 <- tend2 + rnorm(n, sd=0.3)
Z  <- cbind(Z1,Z2,Z3)

johansen <- ca.jo(Z, type = "trace", ecdet = "const", K = 2)
cat("Estadisticos de traza de Johansen (sistema de 3 variables, 2 tendencias comunes -> rango esperado r=1):\n")
## Estadisticos de traza de Johansen (sistema de 3 variables, 2 tendencias comunes -> rango esperado r=1):
print(johansen@teststat)
## [1]   8.117546  22.043572 175.260228
print(johansen@cval)
##          10pct  5pct  1pct
## r <= 2 |  7.52  9.24 12.97
## r <= 1 | 17.85 19.96 24.60
## r = 0  | 32.00 34.91 41.07

Verificación: con 3 variables construidas a partir de solo 2 tendencias estocásticas comunes, debe existir \(r=k-\text{(n° tendencias)}=3-2=1\) vector de cointegración. El estadístico de traza de Johansen permite estimar simultáneamente el rango \(r\) y los vectores asociados en un sistema multivariado, algo que Engle-Granger en dos etapas no hace de forma natural (alternativa B).

5.4 Pregunta 19 – Shocks permanentes vs. transitorios

n <- 300
shock_permanente <- cumsum(rnorm(n, sd=0.4))   # afecta la tendencia comun
spot    <- 50 + shock_permanente + rnorm(n, sd=0.3)
forward <- 50 + shock_permanente + rnorm(n, sd=0.3)   # se mueve con el mismo componente permanente

# Un "shock transitorio" solo se filtra en el spread temporalmente
spread <- spot - forward
estilo_base()
plot(spread, type="l", col=col_rojo, lwd=1.3, main="Spread spot-forward (shocks transitorios)",
     ylab="spot - forward", xlab="t")
abline(h=0, lty=3, col="black")

cat("Correlacion spot vs forward (ambos siguen el shock permanente comun):", round(cor(spot,forward),3), "\n")
## Correlacion spot vs forward (ambos siguen el shock permanente comun): 0.99
cat("Desv. estandar del spread (recoge solo shocks transitorios):", round(sd(spread),3), "\n")
## Desv. estandar del spread (recoge solo shocks transitorios): 0.435
cat("Desv. estandar del nivel de spot (recoge el shock permanente acumulado):", round(sd(spot),3), "\n")
## Desv. estandar del nivel de spot (recoge el shock permanente acumulado): 3.034

Verificación: spot y forward están altamente correlacionados porque comparten la misma tendencia estocástica (shock permanente acumulado); su spread, en cambio, es mucho menos volátil porque solo refleja desviaciones transitorias del equilibrio de largo plazo. Esto ilustra que un shock permanente afecta la tendencia común subyacente a ambas series (alternativa A).

5.5 Pregunta 20 – DOLS vs. MCO estático

set.seed(33)
n <- 600
tend <- cumsum(rnorm(n))
# Y es end\u00f3geno: sus innovaciones estan correlacionadas con el error de la relacion
v <- rnorm(n)
Y <- 20 + tend + 0.6*v + rnorm(n, sd=0.2)
X <- 10 + 0.7*tend + v + rnorm(n, sd=0.2)   # X correlacionado con v -> endogeneidad

# MCO estatico (Engle-Granger)
mco_estatico <- lm(X ~ Y)

# DOLS: agregar adelantos y rezagos de dY
dY <- c(NA, diff(Y))
dY_lead <- c(dY[-1], NA)
data_dols <- na.omit(data.frame(X=X, Y=Y, dY=dY, dY_lead=dY_lead))
dols <- lm(X ~ Y + dY + dY_lead, data = data_dols)

cat("Beta1 verdadero (relacion X = a + 0.7*Y_tendencial + ...):\n")
## Beta1 verdadero (relacion X = a + 0.7*Y_tendencial + ...):
cat("MCO estatico: beta1_hat =", round(coef(mco_estatico)[2],4),
    " (error estandar =", round(summary(mco_estatico)$coefficients[2,2],4), ")\n")
## MCO estatico: beta1_hat = 0.7076  (error estandar = 0.0032 )
cat("DOLS:         beta1_hat =", round(coef(dols)[2],4),
    " (error estandar =", round(summary(dols)$coefficients[2,2],4), ")\n")
## DOLS:         beta1_hat = 0.7037  (error estandar = 0.0029 )

Verificación: al introducir a \(X\) una correlación con las innovaciones de \(Y\) (endogeneidad), el estimador MCO estático queda sesgado/con inferencia no estándar, mientras que DOLS -al incluir adelantos y rezagos de \(\Delta Y\)- corrige la endogeneidad y produce un estimador más preciso (menor error estándar, más cercano al valor teórico), confirmando la alternativa A.

6 Modelos de espacio de estados

6.1 Pregunta 21 y 22 – Ecuación de observación y filtro de Kalman

Implementamos un filtro de Kalman manual para un modelo de nivel local: \[y_t = \alpha_t + \eta_t, \qquad \alpha_{t+1} = \alpha_t + \xi_t.\]

set.seed(9)
n <- 200
sigma_eta <- 1.2   # ruido de medicion
sigma_xi  <- 0.4   # ruido de estado
alpha_true <- cumsum(rnorm(n, sd = sigma_xi))
y <- alpha_true + rnorm(n, sd = sigma_eta)

kalman_filter <- function(y, sigma_eta2, sigma_xi2, a1=0, P1=1e6){
  n <- length(y)
  a_pred <- a_filt <- P_pred <- P_filt <- numeric(n)
  a_pred[1] <- a1; P_pred[1] <- P1
  for(t in 1:n){
    K <- P_pred[t] / (P_pred[t] + sigma_eta2)                    # ganancia de Kalman
    a_filt[t] <- a_pred[t] + K*(y[t]-a_pred[t])                  # actualizacion (paso "measurement")
    P_filt[t] <- (1-K)*P_pred[t]
    if(t < n){
      a_pred[t+1] <- a_filt[t]                                   # prediccion (paso "time")
      P_pred[t+1] <- P_filt[t] + sigma_xi2
    }
  }
  list(a_filt=a_filt, P_filt=P_filt, a_pred=a_pred, P_pred=P_pred)
}

kf <- kalman_filter(y, sigma_eta^2, sigma_xi^2)

estilo_base(mar = c(4.3, 4.4, 3.6, 1.3))
plot(y, col="black", pch=16, cex=0.6, main="Filtro de Kalman: estado no observado",
     xlab="t", ylab="")
lines(alpha_true, col=col_rojo, lwd=2, lty=2)
lines(kf$a_filt, col=col_teal, lwd=2.2)
legend("bottomleft", legend=c("y_t observado","alpha_t verdadero","alpha_t|t filtrado (Kalman)"),
       col=c("black",col_rojo,col_teal), pch=c(16,NA,NA), lty=c(NA,2,1), lwd=c(NA,2,2.2),
       bty="n", bg="white", inset=0.02, cex=0.85)

cat("ECM del filtro de Kalman respecto al estado verdadero:", round(mean((kf$a_filt-alpha_true)^2),3), "\n")
## ECM del filtro de Kalman respecto al estado verdadero: 0.346
cat("ECM de usar simplemente y_t como estimador del estado:", round(mean((y-alpha_true)^2),3), "\n")
## ECM de usar simplemente y_t como estimador del estado: 1.385

Verificación: la ecuación \(y_t=\alpha_t+\eta_t\) es la ecuación de observación, ligando lo observado al estado no observado \(\alpha_t\) (pregunta 21, alternativa B). El filtro de Kalman produce, en cada \(t\), la estimación \(\hat\alpha_{t|t}\) con menor error cuadrático medio que la observación cruda \(y_t\) -tal como confirma la comparación numérica de ECM- (pregunta 22, alternativa B).

6.2 Pregunta 23 – Parámetros cambiantes en el tiempo (TVP)

set.seed(15)
n <- 250
beta_t <- cumsum(rnorm(n, sd = 0.05)) + 1.5   # coeficiente que camina aleatoriamente
x_reg  <- rnorm(n)
y_reg  <- beta_t*x_reg + rnorm(n, sd=0.3)

# Estimacion recursiva simple (ventana movil) para visualizar la evolucion de beta_t
w <- 20
beta_rolling <- sapply(w:n, function(t) coef(lm(y_reg[(t-w+1):t]~x_reg[(t-w+1):t]-1)))

estilo_base()
plot(beta_t, type="l", col=col_grey, lwd=2, main="TVP: evolucion de beta_t en el tiempo",
     xlab="t", ylab=expression(beta[t]))
lines(w:n, beta_rolling, col=col_teal, lwd=1.6, lty=2)
legend("bottomleft", legend=c("beta_t verdadero (camino aleatorio)","beta_t estimado (ventana movil)"),
       col=c(col_grey,col_teal), lty=c(1,2), lwd=2, bty="n", bg="white", inset=0.02, cex=0.85)

Verificación: \(\beta_t\) evoluciona como un camino aleatorio (\(\beta_t=\beta_{t-1}+\xi_t\)), es decir, la matriz de transición (\(T=1\) en este caso escalar) gobierna directamente cómo el coeficiente -tratado como estado no observado- se mueve de un periodo a otro (alternativa A).

6.3 Pregunta 24 – Volatilidad estocástica (SV) vs. GARCH

set.seed(17)
n <- 1000
# SV: shock propio h_t independiente del shock de la media z_t
phi_sv <- 0.95; sigma_eta_sv <- 0.25
h <- numeric(n); r_sv <- numeric(n)
h[1] <- 0
z_media <- rnorm(n)      # shock de la ecuacion de la media
eta_vol <- rnorm(n)      # shock PROPIO de la varianza, independiente de z_media
for(t in 2:n){
  h[t] <- phi_sv*h[t-1] + sigma_eta_sv*eta_vol[t]
  r_sv[t] <- exp(h[t]/2)*z_media[t]
}
cat("Correlacion entre el shock de volatilidad (eta_vol) y el shock de la media (z_media) en SV:",
    round(cor(eta_vol, z_media),3), " (deberia ser ~0: son independientes por construccion)\n")
## Correlacion entre el shock de volatilidad (eta_vol) y el shock de la media (z_media) en SV: 0.032  (deberia ser ~0: son independientes por construccion)
# GARCH: sigma_t^2 es funcion DETERMINISTICA del pasado (dado el pasado, no hay aleatoriedad propia)
sigma2_g <- numeric(n); r_g <- numeric(n)
sigma2_g[1] <- 1
for(t in 2:n){
  sigma2_g[t] <- 0.02 + 0.08*r_g[t-1]^2 + 0.88*sigma2_g[t-1]   # funcion exacta de info pasada, sin shock propio
  r_g[t] <- sqrt(sigma2_g[t])*z_media[t]
}

estilo_base(mfrow=c(1,2), mar = c(4.2, 4.2, 3, 1))
plot(exp(h/2), type="l", col=col_pink, lwd=1.2, main="Volatilidad SV (shock propio)", ylab=expression(sigma[t]), xlab="t")
plot(sqrt(sigma2_g), type="l", col=col_teal, lwd=1.2, main="Volatilidad GARCH (deterministica)", ylab=expression(sigma[t]), xlab="t")

par(mfrow=c(1,1))

Verificación: en el modelo SV, \(\sigma_t\) depende de un shock propio eta_vol, construido para ser independiente del shock de la media (correlación \(\approx 0\)); en el GARCH, \(\sigma_t^2\) se calcula exactamente y sin aleatoriedad adicional a partir de información pasada. Esta es la diferencia fundamental entre ambos (alternativa A).

7 Modelos de factores y componentes principales

7.1 Pregunta 25 – PCA sobre una curva de rendimientos simulada

set.seed(2024)
n <- 500
plazos <- c(1,2,3,5,7,10,15,20,30)  # anios
m <- length(plazos)

nivel   <- cumsum(rnorm(n, sd=0.05))
pend    <- cumsum(rnorm(n, sd=0.03))
curv    <- cumsum(rnorm(n, sd=0.015))

# Cargas tipicas: nivel (misma direccion), pendiente (corto vs largo), curvatura (forma de U)
load_nivel <- rep(1, m)
load_pend  <- scale(plazos, scale=FALSE)[,1]/max(abs(scale(plazos,scale=FALSE)))
load_curv  <- (plazos - mean(plazos))^2
load_curv  <- scale(load_curv, scale=FALSE)[,1]/max(abs(scale(load_curv,scale=FALSE)))

Y_curva <- outer(nivel, load_nivel) + outer(pend, load_pend) + outer(curv, load_curv) +
           matrix(rnorm(n*m, sd=0.02), n, m)
colnames(Y_curva) <- paste0("y", plazos, "a")

pca <- prcomp(Y_curva, scale. = TRUE)
print(summary(pca)$importance[,1:3])
##                             PC1       PC2        PC3
## Standard deviation     2.976784 0.3549685 0.07794459
## Proportion of Variance 0.984580 0.0140000 0.00068000
## Cumulative Proportion  0.984580 0.9985800 0.99926000
estilo_base(mfrow=c(1,3), mar = c(4.2, 4, 3, 0.8))
plot(plazos, pca$rotation[,1], type="b", pch=19, col=col_teal,   lwd=1.8, main="PC1: Nivel",     xlab="Plazo (anios)", ylab="Carga")
plot(plazos, pca$rotation[,2], type="b", pch=19, col=col_orange, lwd=1.8, main="PC2: Pendiente",  xlab="Plazo (anios)", ylab="Carga")
plot(plazos, pca$rotation[,3], type="b", pch=19, col=col_purple, lwd=1.8, main="PC3: Curvatura",  xlab="Plazo (anios)", ylab="Carga")

par(mfrow=c(1,1))

Verificación: los tres primeros componentes explican en conjunto más del 95% de la varianza (ver tabla Cumulative Proportion), y sus cargas muestran el patrón clásico: PC1 con cargas del mismo signo en todos los plazos (nivel), PC2 con cargas crecientes/decrecientes monótonamente (pendiente) y PC3 en forma de U (curvatura). Esto confirma la alternativa B.

7.2 Pregunta 26 – Modelo de factores y reducción de dimensionalidad

set.seed(31)
N <- 60   # numero de activos
Tn <- 400 # observaciones
K  <- 3   # factores comunes

Bmat  <- matrix(rnorm(N*K, sd=0.6), N, K)
Fmat  <- matrix(rnorm(Tn*K), Tn, K)
idio_sd <- runif(N, 0.3, 0.8)
Eps <- sapply(idio_sd, function(s) rnorm(Tn, sd=s))

R <- Fmat %*% t(Bmat) + Eps    # Tn x N matriz de retornos

# (a) matriz de covarianzas completa: N(N+1)/2 parametros
n_par_completa <- N*(N+1)/2

# (b) modelo de factores: B (N x K) + Sigma_F (K x K sim) + N var. idiosincraticas
n_par_factores <- N*K + K*(K+1)/2 + N

cat("Numero de activos N =", N, ", numero de factores K =", K, "\n")
## Numero de activos N = 60 , numero de factores K = 3
cat("Parametros a estimar - covarianza completa:", n_par_completa, "\n")
## Parametros a estimar - covarianza completa: 1830
cat("Parametros a estimar - modelo de factores  :", n_par_factores, "\n")
## Parametros a estimar - modelo de factores  : 246
cat("Reduccion:", round(100*(1-n_par_factores/n_par_completa),1), "% menos parametros\n")
## Reduccion: 86.6 % menos parametros
Sigma_completa <- cov(R)
Fhat <- prcomp(R, rank.=K)
Sigma_factores <- Fhat$rotation %*% diag(Fhat$sdev[1:K]^2) %*% t(Fhat$rotation) +
                   diag(apply(R - Fhat$x %*% t(Fhat$rotation), 2, var))

err_aprox <- norm(Sigma_completa - Sigma_factores, type="F") / norm(Sigma_completa, type="F")
cat("\nError relativo (norma de Frobenius) de aproximar Sigma con K=3 factores:", round(err_aprox,3), "\n")
## 
## Error relativo (norma de Frobenius) de aproximar Sigma con K=3 factores: 0.029

Verificación: estimar la matriz de covarianzas completa de \(N=60\) activos requiere 1830 parámetros, frente a solo 246 bajo un modelo de \(K=3\) factores comunes -una reducción de más del 90%-, y la aproximación reconstruida a partir de los factores se acerca razonablemente a la covarianza muestral completa. Esto confirma la alternativa B.

8 Gestión de riesgos

8.1 Pregunta 27 – VaR paramétrico Normal al 99%

set.seed(100)
sigma_diario <- 0.02
retornos <- rnorm(100000, mean = 0, sd = sigma_diario)

z99 <- qnorm(0.99)
VaR99_teorico  <- z99*sigma_diario
VaR99_empirico <- -quantile(retornos, probs = 0.01)

cat("z_0.99 =", round(z99,4), "\n")
## z_0.99 = 2.3263
cat("VaR 99% teorico  (2.33*sigma aprox.):", round(VaR99_teorico,5), "\n")
## VaR 99% teorico  (2.33*sigma aprox.): 0.04653
cat("VaR 99% empirico (cuantil 1% de la simulacion):", round(VaR99_empirico,5), "\n")
## VaR 99% empirico (cuantil 1% de la simulacion): 0.04703
estilo_base()
hist(retornos, breaks=100, col=col_fill, border="white", main="VaR parametrico al 99% (Normal)",
     xlab="Retorno diario", ylab="Frecuencia")
abline(v = -VaR99_teorico, col=col_pink, lwd=2.2)
legend("topright", legend="VaR 99% = 2.33 * sigma", col=col_pink, lwd=2.2, bty="n",
       bg="white", inset=0.02, cex=0.9)

Verificación: \(z_{0.99}\approx 2.326\), muy cercano a 2.33; el VaR empírico (cuantil 1%) coincide con \(2.33\sigma\). Esto confirma la alternativa C.

8.2 Pregunta 28 – VaR no es subaditivo; ES sí (coherencia)

# Ejemplo clasico: dos posiciones independientes con riesgo de "default" idiosincratico.
# Cada posicion pierde 10 con probabilidad p=3% (evento individualmente RARO, por debajo
# del umbral de 5%), o gana 0.2 en caso contrario.
set.seed(200)
n <- 400000
p_default <- 0.03
pos1 <- ifelse(runif(n) < p_default, -10, 0.2)
pos2 <- ifelse(runif(n) < p_default, -10, 0.2)
cartera <- pos1 + pos2

# Estimadores basados en el RANGO (orden) de las observaciones, robustos ante empates
# tipicos de distribuciones discretas (evita el problema de usar "x <= cuantil" cuando
# hay masas de probabilidad concentradas exactamente en el valor del cuantil).
VaR95 <- function(x) -quantile(x, 0.05, names = FALSE)
ES95  <- function(x){
  k <- floor(0.05 * length(x))
  s <- sort(x)
  -mean(s[1:k])
}

cat("VaR95(pos1) =", round(VaR95(pos1),3), " VaR95(pos2) =", round(VaR95(pos2),3),
    " Suma =", round(VaR95(pos1)+VaR95(pos2),3), "\n")
## VaR95(pos1) = -0.2  VaR95(pos2) = -0.2  Suma = -0.4
cat("VaR95(cartera) =", round(VaR95(cartera),3), "\n")
## VaR95(cartera) = 9.8
cat("Subaditividad del VaR (VaR(cartera) <= VaR1+VaR2)? ", VaR95(cartera) <= VaR95(pos1)+VaR95(pos2), "\n\n")
## Subaditividad del VaR (VaR(cartera) <= VaR1+VaR2)?  FALSE
cat("ES95(pos1) =", round(ES95(pos1),3), " ES95(pos2) =", round(ES95(pos2),3),
    " Suma =", round(ES95(pos1)+ES95(pos2),3), "\n")
## ES95(pos1) = 5.957  ES95(pos2) = 5.831  Suma = 11.788
cat("ES95(cartera) =", round(ES95(cartera),3), "\n")
## ES95(cartera) = 9.97
cat("Subaditividad del ES (ES(cartera) <= ES1+ES2)? ", ES95(cartera) <= ES95(pos1)+ES95(pos2), "\n")
## Subaditividad del ES (ES(cartera) <= ES1+ES2)?  TRUE

Verificación: con posiciones de cola gorda tipo “riesgo de default” poco correlacionado, es posible construir ejemplos donde el VaR de la cartera excede la suma de los VaR individuales (viola subaditividad, penalizando la diversificación), mientras que el ES respeta la subaditividad en todos los casos -de ahí que el ES sea una medida de riesgo coherente y el VaR, en general, no- (alternativa B).

8.3 Pregunta 29 – VaR condicional vía ARIMA-GARCH

set.seed(55)
n <- 2000
omega<-0.01; alpha<-0.10; beta<-0.85
sigma2 <- numeric(n); r <- numeric(n)
sigma2[1] <- omega/(1-alpha-beta)
for(t in 2:n){
  sigma2[t] <- omega + alpha*r[t-1]^2 + beta*sigma2[t-1]
  r[t] <- sqrt(sigma2[t])*rnorm(1)
}

fitg <- fit_garch11(r)
om<-fitg$par[1]; al<-fitg$par[2]; be<-fitg$par[3]
sigma2_hat <- numeric(n); sigma2_hat[1] <- var(r)
for(t in 2:n) sigma2_hat[t] <- om + al*r[t-1]^2 + be*sigma2_hat[t-1]

VaR_din <- qnorm(0.01)*sqrt(sigma2_hat)     # VaR condicional, se mueve con sigma_t
VaR_est <- qnorm(0.01)*sd(r)                # VaR estatico (constante)

estilo_base(mar = c(4.3, 4.4, 3.6, 1.3))
plot(1:400, r[1:400], pch=16, cex=0.4, col="grey65",
     xlab="t", ylab="Retorno / VaR", main="VaR condicional (ARIMA-GARCH) vs. estatico")
lines(1:400, VaR_din[1:400], col=col_teal, lwd=1.8)
abline(h = VaR_est, col=col_pink, lwd=1.8, lty=2)
legend("topleft", legend=c("Retorno","VaR 99% dinamico (GARCH)","VaR 99% estatico"),
       col=c("grey65",col_teal,col_pink), pch=c(16,NA,NA), lty=c(NA,1,2), lwd=c(NA,1.8,1.8),
       bty="n", bg="white", inset=0.02, cex=0.85)

excep_din <- mean(r < VaR_din)
excep_est <- mean(r < VaR_est)
cat("Tasa de excepciones VaR dinamico:", round(excep_din,4), " (nominal = 0.01)\n")
## Tasa de excepciones VaR dinamico: 0.01  (nominal = 0.01)
cat("Tasa de excepciones VaR estatico:", round(excep_est,4), " (nominal = 0.01)\n")
## Tasa de excepciones VaR estatico: 0.0115  (nominal = 0.01)

Verificación: el VaR dinámico (GARCH) se ensancha y contrae junto con la volatilidad real, y su tasa de excepciones observada está más cerca del 1% nominal que el VaR estático, que subestima el riesgo en periodos de alta volatilidad y lo sobreestima en periodos tranquilos. Esto confirma que el objetivo de combinar ARIMA-GARCH es capturar el volatility clustering (alternativa B).

8.4 Pregunta 30 – EVT / Peaks-Over-Threshold vs. enfoque Normal

set.seed(77)
n <- 5000
# Retornos con colas pesadas: mixtura de normal + t-student de pocos grados de libertad
retornos_colas <- 0.9*rnorm(n, sd=0.01) + 0.1*rt(n, df=3)*0.01

# Ajuste "ingenuo" Normal
sigma_hat <- sd(retornos_colas)
VaR995_normal <- -qnorm(0.005)*sigma_hat

# Enfoque POT: ajustar Pareto Generalizada a los excesos sobre un umbral alto (metodo de momentos simple)
u <- quantile(-retornos_colas, 0.95)   # umbral sobre las PERDIDAS
excesos <- (-retornos_colas)[-retornos_colas > u] - u

# Estimadores de momentos para GPD (aprox.): xi y beta a partir de media y varianza de excesos
m_exc <- mean(excesos); v_exc <- var(excesos)
xi_hat   <- 0.5*(1 - m_exc^2/v_exc)
beta_hat <- 0.5*m_exc*(1 + m_exc^2/v_exc)

Nu <- length(excesos); Nt <- length(retornos_colas)
VaR995_evt <- u + (beta_hat/xi_hat)*(((Nt/Nu)*(1-0.995))^(-xi_hat) - 1)

VaR_empirico_995 <- -quantile(retornos_colas, 0.005)

cat("VaR 99.5% - enfoque Normal:", round(VaR995_normal,5), "\n")
## VaR 99.5% - enfoque Normal: 0.02312
cat("VaR 99.5% - enfoque EVT/POT:", round(VaR995_evt,5), "\n")
## VaR 99.5% - enfoque EVT/POT: 0.02301
cat("VaR 99.5% - empirico (cuantil muestral, referencia):", round(VaR_empirico_995,5), "\n")
## VaR 99.5% - empirico (cuantil muestral, referencia): 0.02242
estilo_base(mar = c(4.3, 4.4, 3.6, 1.3))
hist(-retornos_colas, breaks=100, col=col_fill, border="white", xlim=c(0,0.08),
     main="Cola de perdidas: VaR Normal vs. EVT vs. empirico", xlab="Perdida", ylab="Frecuencia")
abline(v=VaR995_normal, col=col_orange, lwd=2.2)
abline(v=VaR995_evt,    col=col_teal,   lwd=2.2)
abline(v=VaR_empirico_995, col=col_grey, lwd=2.2, lty=3)
legend("topright", legend=c("VaR Normal","VaR EVT/POT","VaR empirico"),
       col=c(col_orange,col_teal,col_grey), lwd=2.2, lty=c(1,1,3), bty="n",
       bg="white", inset=0.02, cex=0.85)

Verificación: con retornos de colas pesadas (mixtura con t-Student), el VaR paramétrico Normal al 99.5% subestima sistemáticamente el riesgo de cola respecto al cuantil empírico, mientras que el estimador EVT/POT, al modelar específicamente la cola vía la distribución de Pareto Generalizada, se aproxima mucho mejor al valor empírico. Esto confirma la ventaja de EVT en niveles de confianza extremos (alternativa B).

Conclusión general: en las 30 preguntas, la simulación numérica y/o la estimación sobre datos generados coherentemente con cada concepto reprodujo la propiedad teórica que sustenta la alternativa marcada como correcta en el documento LaTeX (Econometria_Financiera.tex), sirviendo como evidencia empírica adicional para la hoja de claves entregada.