rm(list=ls())
#setwd("D:/亞洲大學上課資料/Portfolio management 2017 Spring")
dat = read.csv("capm.csv",header=T)
# Convert data into time series
library(xts)
## Warning: package 'xts' was built under R version 3.5.3
## Loading required package: zoo
## Warning: package 'zoo' was built under R version 3.5.3
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
dat.xts<-xts(dat[,2:6], order.by= as.Date(dat[,1], "%Y/%m/%d"))
# Sample period
dat.xts.sample<-dat.xts['199311/199811']
# Calculate returns for three stocks and minus daily treasury bill rate
# Daily treasury bill rate can be obtained by converting annual rate into daily rates
riskfree<-dat.xts.sample[,1]/(100*252)
EX_R_sp500 = dat.xts.sample[,2]/lag(dat.xts.sample[,2],1) - 1 - riskfree
EX_R_msft = dat.xts.sample[,3]/lag(dat.xts.sample[,3],1) - 1 - riskfree
EX_R_ge = dat.xts.sample[,4]/lag(dat.xts.sample[,4],1) - 1 - riskfree
EX_R_ford = dat.xts.sample[,5]/lag(dat.xts.sample[,5],1) - 1 - riskfree
# Delete the first row NA data
fit1 = lm(EX_R_msft[-1,]~EX_R_sp500[-1,])
fit2 = lm(EX_R_ge[-1,]~EX_R_sp500[-1,])
fit3 = lm(EX_R_ford[-1,]~EX_R_sp500[-1,])
#options(digits=3)
summary(fit1)
##
## Call:
## lm(formula = EX_R_msft[-1, ] ~ EX_R_sp500[-1, ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.056454 -0.010372 -0.001141 0.009763 0.087598
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.0012503 0.0004745 2.635 0.00852 **
## EX_R_sp500[-1, ] 1.2823959 0.0531064 24.148 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.01692 on 1275 degrees of freedom
## Multiple R-squared: 0.3138, Adjusted R-squared: 0.3133
## F-statistic: 583.1 on 1 and 1275 DF, p-value: < 2.2e-16
summary(fit2)
##
## Call:
## lm(formula = EX_R_ge[-1, ] ~ EX_R_sp500[-1, ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.032201 -0.006673 -0.000135 0.006388 0.036949
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.0003563 0.0002683 1.328 0.184
## EX_R_sp500[-1, ] 1.1982064 0.0300263 39.905 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.009569 on 1275 degrees of freedom
## Multiple R-squared: 0.5554, Adjusted R-squared: 0.555
## F-statistic: 1592 on 1 and 1275 DF, p-value: < 2.2e-16
summary(fit3)
##
## Call:
## lm(formula = EX_R_ford[-1, ] ~ EX_R_sp500[-1, ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.063903 -0.009529 -0.001148 0.008793 0.075596
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.0003010 0.0004421 0.681 0.496
## EX_R_sp500[-1, ] 1.0251141 0.0494769 20.719 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.01577 on 1275 degrees of freedom
## Multiple R-squared: 0.2519, Adjusted R-squared: 0.2513
## F-statistic: 429.3 on 1 and 1275 DF, p-value: < 2.2e-16
# Extract beta for fit1
beta_hat<-rbind(fit1$coefficients, fit2$coefficients, fit3$coefficients)
b_hat<-t(beta_hat)
#degree of freedom
n<-dim(EX_R_msft[-1,])[1]
df = n-2
#residual variance
res_var = c(sum((fit1$residuals)^2), sum((fit2$residuals)^2), sum((fit3$residuals)^2))/df
diagD_hat = diag(res_var)
diagD_hat
## [,1] [,2] [,3]
## [1,] 0.0002864103 0.000000e+00 0.0000000
## [2,] 0.0000000000 9.155882e-05 0.0000000
## [3,] 0.0000000000 0.000000e+00 0.0002486
# covariance matrix by single factor model
Y = cbind(EX_R_msft[-1,], EX_R_ge[-1,], EX_R_ford[-1,])
cov_factor = var(Y)*t(b_hat)%*%b_hat + diag(diagD_hat)
cov_factor
## Close.msft Close.ge Close.ford
## Close.msft 0.0009723001 0.0004517456 0.0004325921
## Close.ge 0.0002568940 0.0003869549 0.0002096391
## Close.ford 0.0003947818 0.0003666803 0.0005975278
#---------------------------------------------------------------------------
# You can also use OLS formula: beta=inv(X'X)X'Y to get the estimated beta
#---------------------------------------------------------------------------
n=dim(EX_R_sp500)[1]
ones = rep(1,n)
X = cbind(ones, EX_R_sp500)
X = as.matrix(X[-1,])
Y = cbind(EX_R_msft, EX_R_ge, EX_R_ford)
Y = as.matrix(Y[-1,])
b_hat.1 = solve(t(X)%*%X)%*%t(X)%*%Y
b_hat.1
## Close.msft Close.ge Close.ford
## ones 0.001250282 0.0003562714 0.0003009948
## Close.sp500 1.282395932 1.1982064408 1.0251141308
b_hat
## [,1] [,2] [,3]
## (Intercept) 0.001250282 0.0003562714 0.0003009948
## EX_R_sp500[-1, ] 1.282395932 1.1982064408 1.0251141308
# follow the formula in the slides
E_hat = Y - X%*%b_hat.1
res_var.1 = diag(t(E_hat)%*%E_hat)/((n-1)-2)
diagD_hat.1 = diag(res_var.1)
diagD_hat.1
## [,1] [,2] [,3]
## [1,] 0.0002864103 0.000000e+00 0.0000000
## [2,] 0.0000000000 9.155882e-05 0.0000000
## [3,] 0.0000000000 0.000000e+00 0.0002486
# covariance matrix by single factor model
cov_factor.1 = var(Y)*t(b_hat)%*%b_hat + diag(diagD_hat)
cov_factor.1
## Close.msft Close.ge Close.ford
## Close.msft 0.0009723001 0.0004517456 0.0004325921
## Close.ge 0.0002568940 0.0003869549 0.0002096391
## Close.ford 0.0003947818 0.0003666803 0.0005975278
cov_factor
## Close.msft Close.ge Close.ford
## Close.msft 0.0009723001 0.0004517456 0.0004325921
## Close.ge 0.0002568940 0.0003869549 0.0002096391
## Close.ford 0.0003947818 0.0003666803 0.0005975278