Examen práctico I
Marián Bravo
2024-09-29
Problema: Predicción de Presión en Función de la Temperatura (Datos Hooker)
El paquete alr4 de R contiene el conjunto de datos Hooker, que incluye dos variables: la temperatura de ebullición medida (en grados Fahrenheit, °F) y la presión del aire medida (en pulgadas de mercurio, inHg), recopiladas por el Dr. Joseph Hooker durante sus expediciones en las montañas del Himalaya. A diferencia de otros estudios, estas mediciones fueron tomadas a mayor altitud, lo que afecta la relación entre la presión y la temperatura de ebullición. El objetivo de este estudio es encontrar un modelo de regresión lineal que permita predecir la presión del aire dada una temperatura de ebullición medida.
Parte 1: Ajuste del Modelo de Regresión
1. Ajuste un modelo de regresión lineal simple que relacione la presión del aire como variable dependiente y la temperatura de ebullición como variable independiente:
Pressure = β0 + β1 × Temperature + e
#install.packages("alr4")
library(alr4)## Warning: package 'alr4' was built under R version 4.2.3## Loading required package: car## Warning: package 'car' was built under R version 4.2.3## Loading required package: carData## Warning: package 'carData' was built under R version 4.2.3## Loading required package: effects## Warning: package 'effects' was built under R version 4.2.3## lattice theme set by effectsTheme()
## See ?effectsTheme for details.data("Hooker")
#?Hooker
# Definiendo la predictora y la respuesta
xs<-Hooker$bp # esta es la predictora (var independiente)
ys<-Hooker$pres # esta es la respuesta (var dependiente)
#Estimamos b1 y b0
SXY<-sum((xs - mean(xs))*(ys - mean(ys)))
SXX<-sum((xs - mean(xs))*(xs - mean(xs)))
(b1 <- SXY/SXX)## [1] 0.4402819(b0<-mean(ys) - (b1*mean(xs)))## [1] -64.41275Graficando el modelo
plot(xs, ys)
abline(a=b0, b=b1, col="red", lty=2)
2. Realice un análisis diagnóstico del modelo ajustado:
- Examine los gráficos de residuos.
##otra manera de hacer el modelo es
ob<-lm(ys ~ xs)
summary(ob)##
## Call:
## lm(formula = ys ~ xs)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.61383 -0.24968 -0.09921 0.26365 0.81232
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -64.412751 1.429165 -45.07 <2e-16 ***
## xs 0.440282 0.007444 59.14 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3563 on 29 degrees of freedom
## Multiple R-squared: 0.9918, Adjusted R-squared: 0.9915
## F-statistic: 3498 on 1 and 29 DF, p-value: < 2.2e-16# 2. Residuales vs valores ajustados
plot(ob$fitted.values, ob$residuals, ylim=c(-4,4))
abline(h=0, lty=2, col="green") 
Es claro que los errores son minimos y a vista de buen cubero podriamos pensar que su medía sí es cero, además más o menos tienen una varianza constante.
- Verifique los supuestos clásicos de la regresión (linealidad, normalidad de los residuos, homocedasticidad).
Linealidad: En la gráfica del punto 1. vemos una clara relación lineal entre los datos, no hay ningún outlier ni datos que afecten de mala manera le modelo.
Normalidad de los residuos:
qqnorm(ob$residuals)
qqline(ob$residuals, col="red", lty=2)
Vemos que la normalidad sí es razonable porque los puntos se acercan a la línea.
Homocedasticidad: Usamos un Scale-location plot
plot(ob, which = 3)
par(mfrow=c(2,2))
plot(ob)
El supuesto de heterocedasticidad sí se cumple porque la línea no es horizontal.
3. Si encuentra evidencia de que el modelo no es adecuado, aplique y compare transformaciones que considere necesarias para mejorar el ajuste y validez del modelo.
Una mejora que se podría hacer al modelo es quitarle algunos puntos que creo que afectan de mala manera al modelo, como el punto 1,2 y 3 que se ven en los graficos de Residual vs.Fitted, en el de Sacle Location y en el de Residuals vs Leverage, el el último hasta se salen de la línea de cook de 0.5.
nxs<-xs[-c(1,2,3)]
nys<-ys[-c(1,2,3)]
nob<-lm(nys ~ nxs)
summary(ob)##
## Call:
## lm(formula = ys ~ xs)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.61383 -0.24968 -0.09921 0.26365 0.81232
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -64.412751 1.429165 -45.07 <2e-16 ***
## xs 0.440282 0.007444 59.14 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3563 on 29 degrees of freedom
## Multiple R-squared: 0.9918, Adjusted R-squared: 0.9915
## F-statistic: 3498 on 1 and 29 DF, p-value: < 2.2e-16summary(nob)##
## Call:
## lm(formula = nys ~ nxs)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.43174 -0.10455 -0.01497 0.13895 0.43426
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -58.624654 1.158218 -50.62 <2e-16 ***
## nxs 0.409444 0.006097 67.16 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2121 on 26 degrees of freedom
## Multiple R-squared: 0.9943, Adjusted R-squared: 0.994
## F-statistic: 4510 on 1 and 26 DF, p-value: < 2.2e-16par(mfrow=c(2,2))
plot(nob)
Vemos que el error estándar es menor en el segundo modelo por lo tanto es más válido, además en el grafico 3 se ve evidencia de que no hay heterocedasticidad y en el gráfico 4 todos los dats caen dentro del cook 0.5. Eso sin considerar que el rcuadrado del nuevo modelo muestra más certeza de que Y tenga relación con X.
Parte 2: Interpretación del Modelo Final
1. Una vez que haya identificado un modelo válido:
Interprete los resultados de la regresión, incluyendo la significancia de los coeficientes, el valor de R2, y otros indicadores relevantes.
Ambos coeficientes tienen un valor p de < 2e-16, lo que indica que son altamente significativos (p < 0.001). Esto sugiere que hay evidencia sólida para rechazar la hipótesis nula de que los coeficientes son iguales a cero. El valor de Adjusted R-squared: 0.994 sugiere un buen ajuste.
Proporcione una interpretación clara de los parámetros del modelo
El parámetro b0 -58.62 nos dice que la presión inicial es -58.62 empieza a crecer conforme crece la temperatura a una proporción de 0.4094.
2. Discuta si el modelo ajustado es adecuado para predecir la presión del aire en función de la temperatura. Si el modelo no es adecuado, ¿por qué? ¿Qué otras mejoras podrían aplicarse?
El modelo parece sí estar bien ajustado aunque algunas mejoras podrían ser aplicarle una transformación logaritmica para ver como se comporta el modelo, aunque d esimple vista no creo que tra tansformación sea completamente necesaria.
#Actualizo b1 y b0
coeficientes <- nob$coefficients
# Asignar b0 y b1
b0 <- coeficientes[1]
b1 <- coeficientes[2] Parte 3: Intervalos de Confianza y Predicción
1. Ahora que tiene un modelo que considera válido, construya intervalos de confianza del 95 % para las predicciones de la presión del aire a las siguientes temperaturas: 185°F, 212°F
# El valor crítico en la distribución t(n-2) es
n<-length(ys)
tcritical95<-qt(0.95, n-2)
hatys<-b0+b1*xs
plot(xs, hatys, xlim=c(180, 220), ylim=c(13,30))
abline(a=b0, b=b1, lty=2, col="red")
points(xs, ys, pch=16, col="skyblue")
xstar<-seq(175, 220, by=1)
hatys_star<-b0 + (b1*xstar)
RSS<-sum((ys - hatys)^2)
hatsig2<-RSS/(n-2)
se_band_pred<-sqrt( hatsig2 * (1 + (1/n) + ((xstar - mean(xs))^2)/SXX) ) #intervalo de confianza para cada punto
# Vamos a crear las bandas de confianza para la predicción
# (observación "solita"):
pred_LO<-hatys_star - tcritical95*se_band_pred
pred_UP<-hatys_star + tcritical95*se_band_pred
# Graficamos las bandas de confianza para la predicción:
lines(xstar, pred_UP, col="blue", lty=2)
lines(xstar, pred_LO, col="blue", lty=2)
Vemos que los unicos puntos que no se ajustan tan bien son los que eliminamos en el nuevo modelo (1,2 y 3)
Los intervalos de confianza al 95% son:
pred1 <- b0+b1*185
pred2 <- b0+b1*212
conf1<-sqrt( hatsig2 * (1 + (1/n) + ((185 - mean(xs))^2)/SXX) )
conf2<-sqrt( hatsig2 * (1 + (1/n) + ((212 - mean(xs))^2)/SXX) )
predlow1<-pred1 - tcritical95*conf1
predup1<-pred1 + tcritical95*conf1
predlow2<-pred2 - tcritical95*conf2
predup2<-pred2 + tcritical95*conf2paste("Con un 95% de confianza podemos decir que la presión del aire cuando la temperatura esta en 185°F está entre (",predlow1, ",", predup1,")")## [1] "Con un 95% de confianza podemos decir que la presión del aire cuando la temperatura esta en 185°F está entre ( 16.3065621083111 , 17.9384800581455 )"paste("Con un 95% de confianza podemos decir que la presión del aire cuando la temperatura esta en 212°F está entre (",predlow2, ",", predup2,")")## [1] "Con un 95% de confianza podemos decir que la presión del aire cuando la temperatura esta en 212°F está entre ( 27.3023690027724 , 29.0526593560284 )"2. Compare los intervalos de confianza con los intervalos de predicción. ¿Qué puede inferir sobre la precisión de las predicciones para estos valores de temperatura?
Los de arriba son los intervalos de confianza de predicción, los intervalos de confianza solos serán más angostos porque estos son itervalos para la media, si nuestros puntos caen adentro de las lineas de confizanda de predicción vamos por buen camino pero un indiador aun más fuerte de que la predicción esta bien es que caiga dentro del intervalo de confianza sobre la media.
Problema: Aplicación del Modelo CAPM usando Regresión Lineal Simple
El Capital Asset Pricing Model (CAPM) es una teoría financiera que describe la relación entre el rendimiento esperado de un activo y su riesgo sistemático en comparación con el mercado en su conjunto. Según el CAPM, el rendimiento esperado de una acción está determinado por su exposición al riesgo del mercado y se puede expresar con la siguiente fórmula:
E(Ri) = Rf + β(E(Rm) − Rf )
El parámetro β es crucial, ya que indica el nivel de riesgo sistemático de una acción en comparación con el mercado. Este parámetro se puede estimar utilizando un modelo de regresión lineal simple, en el cual la variable dependiente es el retorno de la acción y la variable independiente es el retorno del mercado.
Parte 1: Ajuste del Modelo de Regresión CAPM
Elige una acción estadounidense y un índice de mercado relevante (puede ser del mercado o de la industria de la acción seleccionada). ### 1. Descarga los datos históricos de la acción y el índice de mercado para el mismo período de tiempo (2 años y medio), reservando el último medio año para validación. Calcula los retornos diarios usando la diferencia de logaritmos de precios, es decir:
Ri = ln(Pi) − ln(Pi−1) donde Ri es el retorno diario basado en el precio actual Pi y el precio del día anterior Pi−1.
2. Obtén la tasa libre de riesgo para el mismo período. Saqué toda la información de Refinitiv Eikon incluyendo la tasa libre de riesgo que esta en la misma base de datos que los precios.
dt <- read.csv(file = "datacapm.csv", header = TRUE) Datos para el modelo (500 observaciones)
mdt <- head(dt, n = nrow(dt) - 148) ##148 porque dt tiene 648 observaciones y solo quiero 500Retornos de Micorsoft (MSFT.O) y de NASDAQ (NDX):
library(dplyr)## Warning: package 'dplyr' was built under R version 4.2.3##
## Attaching package: 'dplyr'## The following object is masked from 'package:car':
##
## recode## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionmdt$rendMSFT<-log(mdt$MSFT.O)-log(lag(mdt$MSFT.O))
mdt$rendNDX<-log(mdt$.NDX)-log(lag(mdt$.NDX))head(mdt)## Date MSFT.O .NDX USDSOFR. rendMSFT rendNDX
## 1 27/9/2024 428.02 394209.9 4.83 NA NA
## 2 26/9/2024 431.31 394485.8 4.83 0.007657166 0.0006996506
## 3 25/9/2024 432.11 391942.5 4.84 0.001853096 -0.0064681427
## 4 24/9/2024 429.17 385055.1 4.84 -0.006827075 -0.0177284996
## 5 23/9/2024 433.51 385430.4 4.83 0.010061753 0.0009741525
## 6 20/9/2024 435.27 384034.1 4.83 0.004051664 -0.0036294674De los últimos dos años (500 observaciones) porque el último medio año es para validación
3. Ajusta un modelo de regresión lineal simple para estimar el valor de β de la acción.
Nota que el modelo CAPM se puede reescribir como: E(Ri) − Rf = β(E(Rm) − Rf) Esto equivale a realizar una regresión lineal sin intercepto para los retornos en exceso de la acción y el mercado. Calcula Y = Ri−Rf y X = Rm−Rf , y ajusta los siguientes modelos: Y = α + βX + e (SLR) y Y = βX + e (CAPM)
Vemos que β = E(Ri) − Rf / (E(Rm) − Rf)
y <- (mdt$rendMSFT-mdt$USDSOFR.)[-1]
x <- (mdt$rendNDX-mdt$USDSOFR.)[-1]
#les quito el primer dato porque el primer rendimiento es NAMODELO SLR:
slr<-lm(y ~ x)
summary(slr)##
## Call:
## lm(formula = y ~ x)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.059762 -0.006710 -0.000726 0.007172 0.061455
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.002460 0.004309 -0.571 0.568
## x 0.999508 0.000865 1155.517 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.01169 on 497 degrees of freedom
## Multiple R-squared: 0.9996, Adjusted R-squared: 0.9996
## F-statistic: 1.335e+06 on 1 and 497 DF, p-value: < 2.2e-16slrbeta <- slr$coefficients[2]Modelo CAPM:
capm <- lm(y ~ 0 + x)
summary(capm)##
## Call:
## lm(formula = y ~ 0 + x)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.059864 -0.006855 -0.000620 0.007228 0.060466
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## x 0.999998 0.000105 9526 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.01168 on 498 degrees of freedom
## Multiple R-squared: 1, Adjusted R-squared: 1
## F-statistic: 9.075e+07 on 1 and 498 DF, p-value: < 2.2e-16capmbeta <- capm$coefficientspaste("Estimación de β con SLR = ", slrbeta)## [1] "Estimación de β con SLR = 0.999508103261038"paste("Estimación de β con CAPM = ", capmbeta)## [1] "Estimación de β con CAPM = 0.999998282376723"4. Evalúa si el modelo CAPM tiene sentido comparado con el modelo de regresión lineal simple.
Vemos que el p valor de la alpha (b0) del modelo SLR es muy alto lo que indica poca significancia para ese valor lo que no puede llevar a pensar que un modelo de regresión lineal centrado en el origen (CAPM) sería lo mejor. Además el Rajustado en el modelo de CAPM es 1!
5. Interpreta el valor estimado de β. ¿Qué indica en términos del riesgo de la acción en comparación con el mercado? ¿Es la acción más volátil que el mercado (si β > 1) o menos volátil (si β < 1)?
β=1: La acción tiene la misma volatilidad que el mercado. Su rendimiento se moverá en línea con el rendimiento del mercado. β>1: La acción es más volátil que el mercado. Esto significa que si el mercado sube o baja, el rendimiento de la acción tiende a moverse en una magnitud mayor. Por ejemplo, si β=1.5 entonces la acción es 50% más volatil que el mercado. β>1: La acción es menos volátil que el mercado. Por ejemplo, si β=0.8 entonces la acción es 20% mensos volatil que el mercado.
Dado que ya se decidió que la mejor estimación de β es la del modelo CAPM podemos decir que MICROSOFT tiene una volatilidad poco menor pero casi igual a la del mercado. Inclusive se puede ver en la gráfica de Refinitiv de comparación de NASDAQ con MICROSOFT.
Parte 2: Diagnóstico y Análisis del Modelo
Realiza el siguiente análisis una vez ajustado el modelo: ### 1. Verifica los supuestos del modelo de regresión lineal: Realiza un gráfico de residuos para evaluar la presencia de heterocedasticidad. Verifica la normalidad de los residuos utilizando un gráfico Q-Q. Comenta sobre la validez del modelo y su ajuste. ¿Qué aspectos podrían mejorarse?
par(mfrow=c(2,2))
plot(slr)
En las graficas 1, 3 y 4 se ven las distancias salteadas en el eje x porque la tasa libre de riesgo cambia de manera no constante, es decir cualquier día la Reserva Federal puede cambiarla de 4% a 5% y eso es lo que hace cada cierto tiempo.
En el grafico 1 podemos ver que no hay relaicones no lineales en los residuos, en el gráfico 2 vemos que los residuos se distribuyen de manera normal, en la grafica 3 vemos que no hay señales de heterocedasticidad. Una mejora que se puede hacer es quitar algunos outliers que se observan como los puntos 483,233 y 359.
2. Calcula el coeficiente de determinación R2 para los datos de prueba, ¿es muy diferente que el que se obtiene con los datos de entrenamiento?. ¿Qué porcentaje de la variabilidad de los retornos de la acción es explicado por los retornos del mercado?
Los datos de prueba son:
dt$rendMSFT<-log(dt$MSFT.O)-log(lag(dt$MSFT.O))
dt$rendNDX<-log(dt$.NDX)-log(lag(dt$.NDX))
xp <- (dt$rendNDX-dt$USDSOFR.)[c(seq(501, 647, by=1))]yp <- slr$coefficients[1]+slrbeta*xp
ys <- (dt$rendMSFT-dt$USDSOFR.)[c(seq(501, 647, by=1))]
SYY<-sum((ys - mean(ys))^2)
RSS<-sum((ys - yp)^2)
SSreg<-SYY - RSS
(R2<-SSreg/SYY)## [1] 0.9998378El R ajustado que se obtiene con los datos de prueba es practicamente el mismo que se obtiene con los datos de entrenamiento, casi el 100% de la variabilidad de los retornos de la acción es explicado por los retornos del mercado.
3. Al tomar el modelo CAPM, (es decir, una regresión por el origen). ¿El ajuste mejora o empeora?¿Es este cambio significativo?
ypcapm <- capmbeta*xp
RSScapm<-sum((ys - ypcapm)^2)
SSreg<-SYY - RSScapm
(R2capm<-SSreg/SYY)## [1] 0.9998405Vemos que se obtiene un R ajustado ligeramente más grande, lo que sugiere que el modelo CAPM es más acertado y explica mejor el comportamiento de los retornos de la acción.
La diferencia entre los dos Rajustados es de:
(R2capm-R2)## [1] 2.685646e-06Por lo tanto no es demasiada la diferencia entre los dos modelos.
Parte 3: Análisis de Predicción e Intervalos de Confianza
Ahora que has ajustado el modelo, realiza las siguientes tareas: ### 1. Usa el modelo de regresión para predecir el retorno de la acción cuando el retorno del mercado es del 5 %. Calcula un intervalo de predicción del 95% para esta estimación.
2. Grafica los datos y la línea de regresión para el periodo de validación. ¿El modelo proporciona una buena predicción? ¿Puedes incluir las bandas de confianza en el gráfico?
¿Cómo medirías su precisión?
rf <- dt$USDSOFR.[1] ##tasa libre de riesgo de hoy
yp1 <- slr$coefficients[1]+slrbeta*(0.05-rf)
#Como Y = Ri−Rf, entonces Ri =
(re1 <- yp1 + rf) #rendimiento esperado ## (Intercept)
## 0.04989131Para sacar el upper y lower boundary que además nos servirá para graficar
# El valor crítico en la distribución t(n-2) es
n<-length(yp)
tcritical95<-qt(0.95, n-2)
plot(xp, yp, xlim=c(-0.5, 0.1), ylim=c(-1,0.1)) ##observados
abline(a=slr$coefficients[1], b=slrbeta, lty=2, col="red")
points(xp, ys, pch=16, col="skyblue") #predicción
xstar<-seq(-0.5, 0.1, by=0.0005)
hatsig2<-RSS/(n-2)
SXX<-sum((xp - mean(xp))*(xp - mean(xp)))
se_band_pred<-sqrt( hatsig2 * (1 + (1/n) + ((xstar - mean(xp))^2)/SXX) )
hatys_star <- slr$coefficients[1]+slrbeta*xstar
# Vamos a crear las bandas de confianza para la predicción
# (observación "solita"):
pred_LO<-hatys_star - tcritical95*se_band_pred
pred_UP<-hatys_star + tcritical95*se_band_pred
# Graficamos las bandas de confianza para la predicción:
lines(xstar, pred_UP, col="blue", lty=2)
lines(xstar, pred_LO, col="blue", lty=2)
(no se ven todos los datos pero hice un close up para ver mejor que tan ajustados están los datos a las líneas de predicción, vemos que hay pocos que se salen)
Sacando el intervalo de la predicción del rendimiento de la acción:
conf<-sqrt( hatsig2 * (1 + (1/n) + (((0.05-rf) - mean(xp))^2)/SXX) )
predlow<-yp1 - tcritical95*conf
predup<-yp1 + tcritical95*confpaste("Con un 95% de confianza podemos decir que el rendimiento de la acción cuando el rendimiento del mercado esta en 5% está entre (",(predlow+rf)*100, "% ,", (predup+rf)*100,"% )")## [1] "Con un 95% de confianza podemos decir que el rendimiento de la acción cuando el rendimiento del mercado esta en 5% está entre ( 2.99098756708576 % , 6.98727461527104 % )"3. Construye un intervalo de confianza del 95% para β. Interpreta este intervalo en el contexto del modelo CAPM.
#sacamos error estandar del estimador
sehatb1<-sqrt(hatsig2/SXX)
alp<-0.05
t_critic<-qt(alp/2, df=n-2, lower.tail = FALSE)
# Construyendo el IC
(low <- capmbeta - t_critic*sehatb1)## x
## 0.9979083(up <- capmbeta + t_critic*sehatb1)## x
## 1.002088Esto quiere decir que con un 95% de confianza podemos decir que beta esta entre 0.9979083 y 1 (beta por definición no se puede pasar de 1)
Conclusión
Me gustó mucho hacer el modelo de CAPM porque no es siempre en la carrera que vemos ejemplos de la vida real, igual me pareció muy interesante como el modelo centrado en el origen puede ser más certero aunque al pirncipio esto me parecía dificil de entender.