Après le test de différence de moyennes, nous attaquons dans cette session à la régression linéaire simple. Celle-ci permet d’établir le degré d’association entre une variable expliquée (\(y_i\)) et une variable explicative (\(x_i\)). La relation mesurée est matérialisée par une équation linéaire à un facteur (celle d’une droite). On a ainsi : \(y_i=\alpha+\beta.x_i\). Cette approche se distingue du simple coefficient de corrélation dans la mesure où elle implique notamment un sens à la relation. L’évolution de la variable explicative génère l’évolution de la variable expliquée.
L’ensemble est établi au travers des moindres carrés ordinaires (MCO). R propose une fonction simple pour réaliser l’ajustement: lm(). Elle prend comme premier argument l’équation de régression (variable expliquée ~ variable explicative) et comme second (data) le nom de la data frame dans laquelles sont stockées les données (à défaut il faudra indiquer nom_data_frame$ pour le nom de chaque variable dans l’équation). Mais voyons tous ça au travers d’une série d’application.
Tâche n°1 Mesurez au travers d’une regression linéaire simple le degré d’association entre le résultat au test et le ratio enseignés/enseignant
Il s’agit pour nous d’estimer un modèle linéaire sur les données des écoles californiennes pour obtenir approximation de la relation testscr str. Celle-ci prendra la forme:
\[testscr_i=\hat{\alpha}+\hat{\beta}.str_i+\epsilon_i\]
Commençons par charger les packages habituels.
library(tidyverse)
library(haven)
library(gt)Puis, importons les données dans R. Si vous ne les avez pas encore dans votre dossier de travail, je vous renvoie aux deux postes précédents (ici et là).
caschool <- read_dta("caschool.dta")Limitons la base aux seules données utiles str et testscr.
variables <- caschool %>% select(str,testscr)Ceci fait avant d’estimer notre modèle traçons le nuage de points correspondant à l’intersection de nos deux variables. Pour cela, utilisons ggplot() une fonction du tidyerse. Indiquons la data frame dans laquelle les données sont stockées, les variables utilisées en abscisse et un ordonné. Ajoutons avec geom_point() que nous voulons un nuage de points (pour ceux qui voudraient en savoir plus sur ggplot vous trouverez ici une série de posts portant sur son utilisation)
ggplot(data=variables,aes(x=str,y=testscr))+
geom_point()
La forme du nuage nous laisse penser qu’il y aurait une relation négative entre les variables. Il semble qu’en moyenne, plus str est important, moins testscr l’est.
Réalisons le régression pour voir si on a bin un beta négatif.
lm(testscr~str,data=variables)##
## Call:
## lm(formula = testscr ~ str, data = variables)
##
## Coefficients:
## (Intercept) str
## 698.93 -2.28C’est le cas. Voyons ce que cela donne graphiquement en ajoutant à notre précédent graphe une droite reprenant la valeur de paramètres estimées. Pour cela, utilisons le geom_abline().
ggplot(data=variables,aes(x=str,y=testscr))+
geom_point()+
geom_abline(intercept = 698.93,slope = -2.28, color='red')
Voyons si la relation mise au jour apparaît statistiquement significatif. Nous testerons l’hypothèse H0 que le coefficient beta égale à 0 contre H1 beta est différent de 0. Pour cela, le plus simple est d’utiliser la fonction summary() sur notre modèle de régression.
lm(testscr~str,data=variables) %>% summary()##
## Call:
## lm(formula = testscr ~ str, data = variables)
##
## Residuals:
## Min 1Q Median 3Q Max
## -47.727 -14.251 0.483 12.822 48.540
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 698.9330 9.4675 73.825 < 2e-16 ***
## str -2.2798 0.4798 -4.751 2.78e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 18.58 on 418 degrees of freedom
## Multiple R-squared: 0.05124, Adjusted R-squared: 0.04897
## F-statistic: 22.58 on 1 and 418 DF, p-value: 2.783e-06Nous obtenons toute une série de nouvelles informations notamment des tests statistiques sur les paramètres du modèle individuellement (les tests de Student) ou globalement (le test de Fisher sur le dernière ligne). Ici, l’ensemble rejet les hypothèses H0 associées. Si on prend le beta de str, on a pour le test (de Student H0 beta = 0) pratiqué produit une p-value trés trés petite (bien inférieur à 1% 0.01). Cela permet très largement de rejeter l’hypothèse H0. Le coefficient est bien statistiquement différent à 0.
Interprétons rapidement la relation mis en évidence. Tout d’abord, on peut dire si str tend vers 0 alors le score du test tend vers l’ordonnée à l’origine de la droite, le coefficient alpha (la constante du modèle), c’est à dire ici 698.93. Si str = 1, on a \(698.93-2.27/times{1}=696.6\)…
On peut obtenir les prédictions du modèle pour différentes valeurs très facilement en utilisant la fonction predict(). Elle prend comme premier argument le modèle et comme second les valeurs de la variable explicatives pour lesquels vous voulez réaliser la prédiction (dans une data frame). Choisisons 0, 1 , 2 ainsi la moyenne de str et sa médiane.
lm(testscr~str,data=variables) %>%
predict(data.frame(str=c(0,1,2,19.64,19.72))) %>%
data.frame(str=c(0,1,2,19.64,19.72),prediction=.) %>%
`rownames<-`( NULL )## str prediction
## 1 0.00 698.9330
## 2 1.00 696.6531
## 3 2.00 694.3733
## 4 19.64 654.1575
## 5 19.72 653.9751Tâche n°2: Estimez le Beta issu du modèle de marché pour les actions Microsoft, Apple et Facebook (Meta)
La première étape de ce travail consiste à obtenir les données de associé aux cours des actions et de l’indice de marché (ici le S&P500). Pour cela, nous allons utiliser le package quantmod qui propose toute une série de fonction dédiée à des applications de finance de marché.
library(quantmod)La fonction qui nous intéressera ici est getSymbols() qui permet d’obtenir les informations de cours désirées à partir d’un liste de ticker. Ces codes peuvent être retrouvé assez facilement à partir d’une recherche internet. On a ainsi pour Microsoft MSFT, pour Apple AAPL, pour Facebook Meta et pour l’indice S&P500 ^GSPC. Indiquons le site d’extraction yahoo finance ainsi que la période sur laquelle les données doivent courir.
getSymbols(c('MSFT','AAPL','Meta','^GSPC'),src='yahoo',
from="2016-01-01",to="2016-12-31")## [1] "MSFT" "AAPL" "Meta" "^GSPC"L’extraction produit un objet par ticker. Les données sont présentées au format xts qui est spécifique aux séries temporelles. Pour la suite, il est plus simple de les réunir au sein d’une seule data frame classique. Pour chaque ticker, nous établissons une data frame avec une variable date reprenant ce qui était indiqué en nom de ligne dans le xts. Puis, nous fusionnons l’ensemble avec left_join() et suprimons de l’environnement de travail les objets inutils.
t1<-MSFT %>% data.frame %>% data.frame(Date = rownames(.), .) %>%
mutate(date=as.Date(Date))
t2<-AAPL %>% data.frame %>% data.frame(Date = rownames(.), .) %>%
mutate(date=as.Date(Date))
t3<-GSPC %>% data.frame %>% data.frame(Date = rownames(.), .) %>%
mutate(date=as.Date(Date))
t4<-META %>% data.frame %>% data.frame(Date = rownames(.), .) %>%
mutate(date=as.Date(Date))
dat<-left_join(t1,t2) %>% left_join(.,t3) %>% left_join(.,t4)## Joining, by = c("Date", "date")
## Joining, by = c("Date", "date")
## Joining, by = c("Date", "date")rm(t1,t2,t3,t4,MSFT,AAPL,META,GSPC)Limitons l’information aux seuls cours de clôture ajustés (ticker.adjusted).
dat<-dat %>% select(Date,MSFT.Adjusted,AAPL.Adjusted,
META.Adjusted,GSPC.Adjusted)Calculons les rendements journaliers (jours de cotation) simples correspondant aux différentes séries.
\[ rend=\frac{prix_t}{prix_{t-1}}-1\]
dat<-dat %>% mutate(ret_MSFT=MSFT.Adjusted/lag(MSFT.Adjusted)-1,
ret_AAPL=AAPL.Adjusted/lag(AAPL.Adjusted)-1,
ret_META=META.Adjusted/lag(META.Adjusted)-1,
ret_GSPC=GSPC.Adjusted/lag(GSPC.Adjusted)-1)Avant d’aller plus loin sauvegardons des données dans un fichier .csv. Pour cela, il nous suffit d’utiliser la fonction write.table().
write.table(dat,file="return_actions.csv",row.names = FALSE,sep=';')A partir de là, on peut régresser le rendement de chaque titre sur le rendement du marché pour obtenir le Beta du modèle de marché. Pour cela, nous utilisons le fonction lm() et générons les statistiques de présentation du résultat grâce à la fonction summary() associée.
reg1<-summary(lm(ret_MSFT~ret_GSPC,data=dat))
reg1##
## Call:
## lm(formula = ret_MSFT ~ ret_GSPC, data = dat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.071914 -0.004895 0.000246 0.005083 0.047691
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.0001444 0.0006381 0.226 0.821
## ret_GSPC 1.2368714 0.0778015 15.898 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.01009 on 249 degrees of freedom
## (1 observation effacée parce que manquante)
## Multiple R-squared: 0.5037, Adjusted R-squared: 0.5017
## F-statistic: 252.7 on 1 and 249 DF, p-value: < 2.2e-16Le beta correspond au coefficient associé aux rendements du titre. Pour Microsoft, il est de 1,23 et il est statistiquement significative. Le titre amplifie les mouvements de marché.
reg2<-summary(lm(ret_AAPL~ret_GSPC,data=dat))
reg2##
## Call:
## lm(formula = ret_AAPL ~ ret_GSPC, data = dat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.064371 -0.005698 -0.000121 0.005559 0.066085
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.0001047 0.0007666 0.137 0.891
## ret_GSPC 1.0237082 0.0934625 10.953 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.01213 on 249 degrees of freedom
## (1 observation effacée parce que manquante)
## Multiple R-squared: 0.3252, Adjusted R-squared: 0.3224
## F-statistic: 120 on 1 and 249 DF, p-value: < 2.2e-16Pour Apple, il est de 1,02 et il est statistiquement significative. Le titre suit simplement les mouvements de marché.
reg3<-summary(lm(ret_META~ret_GSPC,data=dat))
reg3##
## Call:
## lm(formula = ret_META ~ ret_GSPC, data = dat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.051670 -0.006175 -0.000512 0.005760 0.149045
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.0001248 0.0009768 0.128 0.898
## ret_GSPC 1.0933911 0.1190926 9.181 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.01545 on 249 degrees of freedom
## (1 observation effacée parce que manquante)
## Multiple R-squared: 0.2529, Adjusted R-squared: 0.2499
## F-statistic: 84.29 on 1 and 249 DF, p-value: < 2.2e-16Pour Facebook, il est de 1,09 et il est statistiquement significative. Le titre suit les mouvements de marché.
Réunissons le tout dans un tableau unique. Etablissons une data frame incluant les betas, les t stat associé et les p-values.
tab<-data.frame(c('Microsoft','Apple','Facebook'),
rbind(round(reg1$coefficients[2,c(1,3,4)],digits=4),
round(reg2$coefficients[2,c(1,3,4)],digits=4),
round(reg3$coefficients[2,c(1,3,4)],digits=4)))
colnames(tab)<-c(' ','Beta','t stat', 'p value')
tab$`p value`<-format(tab$`p value`,nsmall = 4)
tab## Beta t stat p value
## 1 Microsoft 1.2369 15.8978 0.0000
## 2 Apple 1.0237 10.9531 0.0000
## 3 Facebook 1.0934 9.1810 0.0000Mettons le résultat sous la forme d’un tableau (gt()) et mettons l’ensemble en page.
tab_<-gt(tab) %>%
tab_header(title='Beta issu du modèle de marché',
subtitle = "Estimé sur 2016 en référence à l'indice S&P500")
tab_| Beta issu du modèle de marché | |||
| Estimé sur 2016 en référence à l'indice S&P500 | |||
| Beta | t stat | p value | |
|---|---|---|---|
| Microsoft | 1.2369 | 15.8978 | 0.0000 |
| Apple | 1.0237 | 10.9531 | 0.0000 |
| 1.0934 | 9.1810 | 0.0000 | |
Importons le résultat dans Word via gtsave().
gtsave(tab_,'Betas_v1.docx')