#Importing data
library(quantmod)## Loading required package: xts## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric## Loading required package: TTR## Version 0.4-0 included new data defaults. See ?getSymbols.library(tidyquant)## Loading required package: lubridate##
## Attaching package: 'lubridate'## The following object is masked from 'package:base':
##
## date## Loading required package: PerformanceAnalytics##
## Attaching package: 'PerformanceAnalytics'## The following object is masked from 'package:graphics':
##
## legend## Loading required package: tidyverse## -- Attaching packages --------------------------------------- tidyverse 1.2.1 --## v ggplot2 3.1.0 v purrr 0.3.2
## v tibble 2.1.1 v dplyr 0.8.0.1
## v tidyr 0.8.3 v stringr 1.4.0
## v readr 1.3.1 v forcats 0.4.0## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x lubridate::as.difftime() masks base::as.difftime()
## x lubridate::date() masks base::date()
## x dplyr::filter() masks stats::filter()
## x dplyr::first() masks xts::first()
## x lubridate::intersect() masks base::intersect()
## x dplyr::lag() masks stats::lag()
## x dplyr::last() masks xts::last()
## x lubridate::setdiff() masks base::setdiff()
## x lubridate::union() masks base::union()library(lubridate)
library(timetk)
library(plyr)## -------------------------------------------------------------------------## You have loaded plyr after dplyr - this is likely to cause problems.
## If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
## library(plyr); library(dplyr)## -------------------------------------------------------------------------##
## Attaching package: 'plyr'## The following objects are masked from 'package:dplyr':
##
## arrange, count, desc, failwith, id, mutate, rename, summarise,
## summarize## The following object is masked from 'package:purrr':
##
## compact## The following object is masked from 'package:lubridate':
##
## herelibrary(quadprog)
library(lpSolve)
tickers<-c("SPY", "QQQ", "EEM", "IWM", "EFA", "TLT", "IYR", "GLD")
data.env<-new.env()
l_ply(tickers, function(sym) try(getSymbols(sym, from="2010-01-01", env=data.env), silent = T))## 'getSymbols' currently uses auto.assign=TRUE by default, but will
## use auto.assign=FALSE in 0.5-0. You will still be able to use
## 'loadSymbols' to automatically load data. getOption("getSymbols.env")
## and getOption("getSymbols.auto.assign") will still be checked for
## alternate defaults.
##
## This message is shown once per session and may be disabled by setting
## options("getSymbols.warning4.0"=FALSE). See ?getSymbols for details.stocks<- tickers[tickers %in% ls(data.env)]
data<- xts()
for(i in seq_along(stocks)){
symbol<-stocks[i]
data <- merge(data, Ad(get(symbol, envir=data.env)))
}
data<- xts(coredata(data), order.by = as.Date(index(data)))
head(data)## SPY.Adjusted QQQ.Adjusted EEM.Adjusted IWM.Adjusted
## 2010-01-04 94.13087 42.07555 35.67281 56.09659
## 2010-01-05 94.38007 42.07555 35.93174 55.90370
## 2010-01-06 94.44649 41.82175 36.00691 55.85109
## 2010-01-07 94.84521 41.84896 35.79811 56.26318
## 2010-01-08 95.16080 42.19340 36.08208 56.57005
## 2010-01-11 95.29370 42.02117 36.00691 56.34210
## EFA.Adjusted TLT.Adjusted IYR.Adjusted GLD.Adjusted
## 2010-01-04 43.40944 68.21629 32.26275 109.80
## 2010-01-05 43.44770 68.65686 32.34021 109.70
## 2010-01-06 43.63134 67.73775 32.32613 111.51
## 2010-01-07 43.46300 67.85172 32.61488 110.82
## 2010-01-08 43.80734 67.82134 32.39655 111.37
## 2010-01-11 44.16698 67.44915 32.55150 112.85#*****************************************************************
# 2. Calculate weekly and monthly returns using log returns
#******************************************************************
data.weekly<-data[endpoints(data, on="weeks", k= 1),]
lrets <- diff(log(data.weekly), lag= 1)
logretw <- na.omit(lrets)
head(logretw)## SPY.Adjusted QQQ.Adjusted EEM.Adjusted IWM.Adjusted
## 2010-01-15 -0.008150151 -0.01515199 -0.02936205 -0.01310461
## 2010-01-22 -0.039762846 -0.03755592 -0.05739643 -0.03110044
## 2010-01-29 -0.016805790 -0.03151480 -0.03415430 -0.02659380
## 2010-02-05 -0.006820688 0.00443036 -0.02861874 -0.01407296
## 2010-02-12 0.012855091 0.01798540 0.03278968 0.02909818
## 2010-02-19 0.028289105 0.02415713 0.02415957 0.03288520
## EFA.Adjusted TLT.Adjusted IYR.Adjusted GLD.Adjusted
## 2010-01-15 -0.003499458 1.984899e-02 -0.006324176 -0.004589865
## 2010-01-22 -0.057354334 1.005025e-02 -0.042683035 -0.033851813
## 2010-01-29 -0.026141733 3.363481e-03 -0.008483892 -0.011354685
## 2010-02-05 -0.019238669 -5.439768e-05 0.003218778 -0.012153577
## 2010-02-12 0.005230991 -1.965291e-02 -0.007602877 0.022294528
## 2010-02-19 0.022734752 -8.238338e-03 0.048966431 0.022447943data.monthly<-data[endpoints(data, on="months", k=1),]
lrets <- diff(log(data.monthly), lag=1)
logretm <- na.omit(lrets)
head(logretm)## SPY.Adjusted QQQ.Adjusted EEM.Adjusted IWM.Adjusted
## 2010-02-26 0.03071820 0.04501028 0.017607768 0.04377887
## 2010-03-31 0.05909806 0.07428098 0.077987148 0.07909460
## 2010-04-30 0.01535186 0.02217746 -0.001663121 0.05523082
## 2010-05-28 -0.08278913 -0.07679835 -0.098645456 -0.07835791
## 2010-06-30 -0.05312741 -0.06161692 -0.014084732 -0.08059607
## 2010-07-30 0.06606925 0.07006902 0.103751388 0.06514060
## EFA.Adjusted TLT.Adjusted IYR.Adjusted GLD.Adjusted
## 2010-02-26 0.00266388 -0.003430628 0.05313394 0.032223420
## 2010-03-31 0.06189853 -0.020788116 0.09302108 -0.004396042
## 2010-04-30 -0.02844674 0.032678280 0.06192350 0.057168648
## 2010-05-28 -0.11870236 0.049822170 -0.05851430 0.030056874
## 2010-06-30 -0.02083487 0.056360158 -0.04782665 0.023280092
## 2010-07-30 0.10984393 -0.009508924 0.08988416 -0.052210719#*****************************************************************
# 3. Convert monthly returns into tibble format.
#******************************************************************
tibble::as_tibble(logretm, rownames = 'date')## Warning: Calling `as_tibble()` on a vector is discouraged, because the behavior is likely to change in the future. Use `tibble::enframe(name = NULL)` instead.
## This warning is displayed once per session.## # A tibble: 113 x 9
## date SPY.Adjusted QQQ.Adjusted EEM.Adjusted IWM.Adjusted EFA.Adjusted
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 2010~ 0.0307 0.0450 0.0176 0.0438 0.00266
## 2 2010~ 0.0591 0.0743 0.0780 0.0791 0.0619
## 3 2010~ 0.0154 0.0222 -0.00166 0.0552 -0.0284
## 4 2010~ -0.0828 -0.0768 -0.0986 -0.0784 -0.119
## 5 2010~ -0.0531 -0.0616 -0.0141 -0.0806 -0.0208
## 6 2010~ 0.0661 0.0701 0.104 0.0651 0.110
## 7 2010~ -0.0460 -0.0527 -0.0329 -0.0774 -0.0387
## 8 2010~ 0.0858 0.124 0.111 0.117 0.0951
## 9 2010~ 0.0375 0.0615 0.0297 0.0406 0.0373
## 10 2010~ 0 -0.00173 -0.0295 0.0343 -0.0494
## # ... with 103 more rows, and 3 more variables: TLT.Adjusted <dbl>,
## # IYR.Adjusted <dbl>, GLD.Adjusted <dbl>#*****************************************************************
# 4. Download Fama French 3 factors data and change
# to digit numbers (not in percetage):
#******************************************************************
library(tidyverse)
library(dplyr)
retdata = read_csv('FFdata3.csv')## Parsed with column specification:
## cols(
## Date = col_double(),
## `Mkt-RF` = col_double(),
## SMB = col_double(),
## HML = col_double(),
## RF = col_double()
## )## Warning: 8 parsing failures.
## row col expected actual file
## 1115 Date a double Annual Factors: January-December 'FFdata3.csv'
## 1115 NA 5 columns 1 columns 'FFdata3.csv'
## 1116 Mkt-RF a double Mkt-RF 'FFdata3.csv'
## 1116 SMB a double SMB 'FFdata3.csv'
## 1116 HML a double HML 'FFdata3.csv'
## .... ...... ......... ................................ .............
## See problems(...) for more details.glimpse(retdata)## Observations: 1,209
## Variables: 5
## $ Date <dbl> 192607, 192608, 192609, 192610, 192611, 192612, 19270...
## $ `Mkt-RF` <dbl> 2.96, 2.64, 0.36, -3.24, 2.53, 2.62, -0.06, 4.18, 0.1...
## $ SMB <dbl> -2.30, -1.40, -1.32, 0.04, -0.20, -0.04, -0.56, -0.10...
## $ HML <dbl> -2.87, 4.19, 0.01, 0.51, -0.35, -0.02, 4.83, 3.17, -2...
## $ RF <dbl> 0.22, 0.25, 0.23, 0.32, 0.31, 0.28, 0.25, 0.26, 0.30,...colnames(retdata)[2]<- 'Mkt_RF'
data <- read.table("FFdata3.csv", sep="\t", dec=",", row.names=1)
stock.rets<-retdata %>% select(c(2, 3, 4))/100
glimpse(stock.rets)## Observations: 1,209
## Variables: 3
## $ Mkt_RF <dbl> 0.0296, 0.0264, 0.0036, -0.0324, 0.0253, 0.0262, -0.000...
## $ SMB <dbl> -0.0230, -0.0140, -0.0132, 0.0004, -0.0020, -0.0004, -0...
## $ HML <dbl> -0.0287, 0.0419, 0.0001, 0.0051, -0.0035, -0.0002, 0.04...N <- dim(stock.rets)[1]
#Mkt.RF<-retdata %>% select(2)/100
fit = lm(formula = cbind(Mkt_RF, SMB, HML)~Mkt_RF, data = stock.rets)
bbeta = as.matrix(fit$coefficients)
bbeta = as.matrix(bbeta[-1,])
bbeta## [,1]
## Mkt_RF 1.0000000
## SMB 0.2457757
## HML 0.1481739sigeps = crossprod(fit$residuals)/(N-2)
# sigeps = as.matrix(var(fit$residuals)) # you can use this way too
sigeps = diag(diag(sigeps))
sigeps## [,1] [,2] [,3]
## [1,] 1.420708e-34 0.000000000 0.000000000
## [2,] 0.000000e+00 0.002045851 0.000000000
## [3,] 0.000000e+00 0.000000000 0.002632401Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.