Tâche n°1 Reprenez les données sur les actions, estimez le modèle de marché pour chacun d’eux, réalisez un interval de confiance à 95% sur le beta, testez sa valeur par rapport à 1.

Comme à chaque fois, nous allons commencer par charger les packages qui vont être utilisés. On y retrouve les habituels tidyverse, haven, gt, mais aussi car, un nouveau package qui inclue des fonctions permettant la réalisation de tests statistiques complémentaires aux régressions.

library(tidyverse)
library(haven)
library(gt)
library(car)

On peut alors importer les données concernant les rendements des actions et de l’indice de marché que nous avons créé dans le post précédent. Si nécessaire vous pouvez le télécharger ici.

return_actions <- read_delim("return_actions.csv",delim = ";")

## Rows: 252 Columns: 9
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ";"
## dbl  (8): MSFT.Adjusted, AAPL.Adjusted, META.Adjusted, GSPC.Adjusted, ret_MS...
## date (1): Date
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

Ceux-ci fait. On peut ré-estimer le modèle de marché pour chaque titre (régression du rendement de l’action sur l’indice de marché en utilisant la fonction lm()) et stocker l’information dans des objets dédiés (portant ici le noms des actions).

Ms<-lm(ret_MSFT~ret_GSPC,data=return_actions)
Ap<-lm(ret_AAPL~ret_GSPC,data=return_actions)
Fb<-lm(ret_META~ret_GSPC,data=return_actions)

Pour obtenir l’interval de confiance au niveau désiré (95%), il suffit de mobiliser la fonction confint(). Elle prend comme premier argument le modèle de régession et comme second (level) le niveau de confiance retenu. Produisons-la pour Microsoft (Ms).

confint(Ms,level=0.95) 

##                    2.5 %      97.5 %
## (Intercept) -0.001112434 0.001401162
## ret_GSPC     1.083637524 1.390103522

On peut lire que le Beta de Microsoft à 95% de chance de se trouver entre 1.083 et 1.390. Néanmoins, l’ensemble reste peut lisible du fait de l’inclusion de l’intervalle de confiance de la constante et du nombre de chiffres après la virgule affiché. Réduisons tout cela.

confint(Ms,parm=2,level=0.95) %>% round(digits = 4)

##           2.5 % 97.5 %
## ret_GSPC 1.0836 1.3901

Créons un tableau incluant l’ensemble des intervalles de confiance des beta des titre.

tb<-rbind(confint(Ms,parm=2,level=0.95),
          confint(Ap,parm=2,level=0.95),
          confint(Fb,parm=2,level=0.95)) %>% round(digits=4)
rownames(tb)<-NULL
tb

##       2.5 % 97.5 %
## [1,] 1.0836 1.3901
## [2,] 0.8396 1.2078
## [3,] 0.8588 1.3279

Incluons-les dans le tableau des beta tel que nous l’avions réalisé dans le post précédent.

reg1<-summary(Ms)
reg2<-summary(Ap)
reg3<-summary(Fb)
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)),
                tb) 
colnames(tab)<-c(' ','Beta','t stat', 'p value','2,5%','97,5%')
tab$`p value`<-format(tab$`p value`,nsmall = 4)
tab

##               Beta  t stat p value   2,5%  97,5%
## 1 Microsoft 1.2369 15.8978  0.0000 1.0836 1.3901
## 2     Apple 1.0237 10.9531  0.0000 0.8396 1.2078
## 3  Facebook 1.0934  9.1810  0.0000 0.8588 1.3279

Ajoutons un titre, un sous-titre ainsi qu’un chapeau sur les colonnes correspondant à l’intervalle de confiance.

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_spanner('Interval de confiance à 95%',columns = c('2,5%','97,5%'))
tab_

Transférons le résultat dans un document Word.

gtsave(tab_,'Betas_conf_v1.docx')

Il nous reste à réaliser les test des beta contre la valeur 1. Les tests dans les modèles de régression sont réalisés contre 0. Pour cela, on peut utiliser la fonction linearHypothesis() du package car. Elle permet de mettre en oeuvre un test de restriction spécification de Fisher. Ses arguments sont le modèle de régression et la restriction (ici ret_GSPC=1).

linearHypothesis(Ms, "ret_GSPC = 1")

## Linear hypothesis test
## 
## Hypothesis:
## ret_GSPC = 1
## 
## Model 1: restricted model
## Model 2: ret_MSFT ~ ret_GSPC
## 
##   Res.Df      RSS Df  Sum of Sq      F  Pr(>F)   
## 1    250 0.026314                                
## 2    249 0.025370  1 0.00094442 9.2693 0.00258 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pratiqué pour Microsoft, on rejette l’hypothèse H1 d’égalité du coefficient à 1 (que le modèle restreint est meilleur que le modèle de base). Appliquons la même procédure qui a un beta plus proche de 1.

linearHypothesis(Ap, "ret_GSPC = 1")

## Linear hypothesis test
## 
## Hypothesis:
## ret_GSPC = 1
## 
## Model 1: restricted model
## Model 2: ret_AAPL ~ ret_GSPC
## 
##   Res.Df      RSS Df  Sum of Sq      F Pr(>F)
## 1    250 0.036621                            
## 2    249 0.036612  1 9.4582e-06 0.0643    0.8

Cette fois ci, on ne peut rejeté H0. Le beta de Apple n’est pas statistiquement différent de 1.

La procédure rempli ses objectifs. Néanmoins, on pourrait préférer utiliser un test de Student comme pour l’hypothèse classique de beta égale 0 que l’on retrouve dans les résumés des modèles de régression. Pour cela, le plus simple (pour moi) est de créer soit même sa fonction de test. Rappelons la forme de statistique de Sudent:

\[Beta=\frac{\hat{\beta}-1}{\sigma/\sqrt{n}}\] où

\[\text{Standard error}=\sigma/\sqrt{n}\]

La fonction pt() permet de calculer la fonction de densité cumulative (la fonction de répartition) de la loi de Student. Mobilisons-la pour calculer la p value du test.

ttest <- function(reg, coefnum, val){
  co <- coef(summary(reg))
  tstat <- (co[coefnum,1]-val)/co[coefnum,2]
  pval<- 2 * pt(abs(tstat), reg$df.residual, lower.tail = FALSE)
  ret<-round(c(tstat,pval),digits=4)
  names(ret)<-c('tstat','pval')
  return(ret)}

Appliquons notre fonction à la régression sur Microsoft (Ms) d’abord en testant beta égal à 0 pour comparer le résultat avec celui obtenu avec la fonction summary(), puis en testant beta égal à 1.

ttest(Ms,2,0)

##   tstat    pval 
## 15.8978  0.0000

ttest(Ms,2,1)

##  tstat   pval 
## 3.0445 0.0026

Notre fonction maison produit bien le résultat attendu. Appliquons-la à toute nos actions. Agrégeons l’ensemble pour l’ajouter à notre tableau des beta.

tc<-rbind(ttest(Ms,2,1),ttest(Ap,2,1),ttest(Fb,2,1))
tc

##       tstat   pval
## [1,] 3.0445 0.0026
## [2,] 0.2536 0.8000
## [3,] 0.7842 0.4337

Créons la data frame correspondant.

tab<-data.frame(tab,tc) 
colnames(tab)<-c(' ','Beta','t stat', 'p value','2,5%','97,5%',
                 't_stat', 'p_value')
tab

##               Beta  t stat p value   2,5%  97,5% t_stat p_value
## 1 Microsoft 1.2369 15.8978  0.0000 1.0836 1.3901 3.0445  0.0026
## 2     Apple 1.0237 10.9531  0.0000 0.8396 1.2078 0.2536  0.8000
## 3  Facebook 1.0934  9.1810  0.0000 0.8588 1.3279 0.7842  0.4337

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_spanner('Beta=0',columns = c('t stat','p value')) %>% 
  tab_spanner('Interval de confiance à 95%',columns = c('2,5%','97,5%')) %>% 
  tab_spanner('Beta=1',columns = c('t_stat','p_value')) 
tab_

Transférons le tableau produit dans Word.

gtsave(tab_,'Betas_conf_0_1_v1.docx')

Tâche n°2 Revenez sur les cas des écoles californiennes, incluez dans la régresion les variables d’éléves apprenant l’anglais (el_pct), étudiez les liens entre les variables

Commençons par recharger les données associés aux écoles (le fichier caschool). Si vous ne les avez pas encore vous pouvez les trouver en lien ici.

caschool <- read_dta("caschool.dta")

Une fois les données importées, créons deux modèles de régression le premier expliquant testscr uniquement avec str (reg1) et le secons l’expliquant à la fois par str et el_pct.

reg1<-lm(testscr~str,data=caschool)
reg2<-lm(testscr~str+el_pct,data=caschool)

Inspectons le résultat de la première régression via summary().

summary(reg1)

## 
## Call:
## lm(formula = testscr ~ str, data = caschool)
## 
## 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-06

Le coefficient associé à str est de -2.2798 et il est statistiquement différent de zéro. Voyons ce qu’il en est dans le second modèle.

summary(reg2)

## 
## Call:
## lm(formula = testscr ~ str + el_pct, data = caschool)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -48.845 -10.240  -0.308   9.815  43.461 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 686.03225    7.41131  92.566  < 2e-16 ***
## str          -1.10130    0.38028  -2.896  0.00398 ** 
## el_pct       -0.64978    0.03934 -16.516  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 14.46 on 417 degrees of freedom
## Multiple R-squared:  0.4264, Adjusted R-squared:  0.4237 
## F-statistic:   155 on 2 and 417 DF,  p-value: < 2.2e-16

La valeur du coefficient de str a changé. Elle est passée à -1.1013 mais est toujours statistiquement significatif. Le coefficient de el_pct est de -0.6497 et est statistiquement significatif.

Pour bien comprendre comment l’inclusion d’el_pct permet d’ajuster l’effet de str sur testscr, régressons str par el_pct. Cela nous indique l’effet de el_pct sur str. Sortons les résidus de la régression autrement dit la part de str non expliquée par el_pct.

rega<-lm(str~el_pct,data=caschool)
res_a<-rega$residuals
head(res_a)

##          1          2          3          4          5          6 
## -1.4444089  2.1013692 -1.2194820 -1.9771752 -0.9320077  1.8310406

Régressons de la même manière testscr par el_pct et sortons le résidus, la part de testscr non expliqué par el_pct.

regb<-lm(testscr~el_pct,data=caschool)
res_b<-regb$residuals
head(res_b)

##           1           2           3           4           5           6 
##  26.0605447  -0.4632982  -1.0047804 -17.0394308 -14.5888014 -50.8611787

Régressons ce dernier résidu par le premier en omettant la constante. Cela se fait en ajoutant -1 aux variables explicatives.

lm(res_b~res_a-1)

## 
## Call:
## lm(formula = res_b ~ res_a - 1)
## 
## Coefficients:
##  res_a  
## -1.101

On retrouve retrouve le même coefficient que dans notre modèle incluant les deux variables pour expliquer testscr. L’inclusion de la seconde variable permet bien de corriger l’effet de str sur testscr de la partie de cet effet qui passe par el_pct.