firm_data1 = read.csv('data.of3firm 2.csv')
str(firm_data1)
## 'data.frame': 12 obs. of 4 variables:
## $ Date : Factor w/ 12 levels "2018/10/1","2018/11/1",..: 4 5 6 7 8 9 1 2 3 10 ...
## $ EPS.BA : num 327 345 331 351 338 ...
## $ EPS.AIR: num 95.8 97.3 100.3 106 106.3 ...
## $ EPS.LMT: num 312 306 289 319 314 ...
firm_data1$date
## NULL
library(xts)
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
library(PerformanceAnalytics)
##
## Attaching package: 'PerformanceAnalytics'
## The following object is masked from 'package:graphics':
##
## legend
date1 = as.Date(firm_data1[,1], "%Y/%m/%d")
#convert firm_data1 into time series data: xts
firm_data1.xts = as.xts(firm_data1[,-1], order.by = date1)
firm.data1<-coredata(firm_data1.xts)
summary(firm.data1)
## EPS.BA EPS.AIR EPS.LMT
## Min. :320.9 Min. : 83.96 Min. :260.0
## 1st Qu.:336.0 1st Qu.: 96.93 1st Qu.:289.6
## Median :348.3 Median :100.24 Median :303.2
## Mean :356.6 Mean :101.82 Mean :301.9
## 3rd Qu.:371.6 3rd Qu.:106.74 3rd Qu.:312.7
## Max. :437.8 Max. :117.90 Max. :341.0
skewness(firm.data1)
## EPS.BA EPS.AIR EPS.LMT
## Skewness 1.309743 -0.02995509 -0.1462904
rbind(apply(firm.data1, 2, summary),
apply(firm.data1, 2, skewness),
apply(firm.data1, 2, kurtosis))
## EPS.BA EPS.AIR EPS.LMT
## Min. 320.887939 83.95999900 259.9671630
## 1st Qu. 336.045677 96.93331950 289.5628125
## Median 348.252762 100.24000150 303.1663055
## Mean 356.593791 101.81860733 301.8978118
## 3rd Qu. 371.610069 106.74000150 312.6891483
## Max. 437.760773 117.90000200 340.9867550
## 1.309743 -0.02995509 -0.1462904
## 1.222041 -0.26193231 0.4379288
#====================================================================================
# IF you know the ticker of stocks, then you can
# download directly from yahoo finance
#=====================================================================================
library(plyr)
library(quantmod)
## Loading required package: TTR
## Version 0.4-0 included new data defaults. See ?getSymbols.
tickers<-c("AIR.PA", "BA", "LMT")
data.env<-new.env()
# here we use l_ply so that we don't double save the data
# getSymbols() does this already so we just want to be memory efficient
# go through every stock and try to use getSymbols()
l_ply(tickers, function(sym) try(getSymbols(sym, 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.
# now we only want the stocks that got stored from getSymbols()
# basically we drop all "bad" tickers
stocks <- tickers[tickers %in% ls(data.env)]
# now we just loop through and merge our good stocks
# if you prefer to use an lapply version here, that is also fine
# since now we are just collecting all the good stock xts() objects
data <- xts()
# i=1
for(i in seq_along(stocks)) {
symbol <- stocks[i]
data <- merge(data, Ad(get(symbol, envir=data.env)))
}
head(data)
## AIR.PA.Adjusted BA.Adjusted LMT.Adjusted
## 2007-01-02 08:00:00 21.50374 NA NA
## 2007-01-03 08:00:00 21.34833 65.44486 62.40288
## 2007-01-04 08:00:00 20.95572 65.70907 62.23996
## 2007-01-05 08:00:00 21.02116 65.43017 62.46399
## 2007-01-08 08:00:00 20.72670 65.27605 63.59760
## 2007-01-09 08:00:00 20.89846 64.58615 63.53651
str(data)
## An 'xts' object on 2007-01-02 08:00:00/2019-04-01 08:00:00 containing:
## Data: num [1:3162, 1:3] 21.5 21.3 21 21 20.7 ...
## - attr(*, "dimnames")=List of 2
## ..$ : NULL
## ..$ : chr [1:3] "AIR.PA.Adjusted" "BA.Adjusted" "LMT.Adjusted"
## Indexed by objects of class: [POSIXct,POSIXt] TZ:
## xts Attributes:
## NULL
# convert POSIXct into date series
data<-xts(coredata(data), order.by = as.Date(index(data), tz=""))
head(data)
## AIR.PA.Adjusted BA.Adjusted LMT.Adjusted
## 2007-01-02 21.50374 NA NA
## 2007-01-03 21.34833 65.44486 62.40288
## 2007-01-04 20.95572 65.70907 62.23996
## 2007-01-05 21.02116 65.43017 62.46399
## 2007-01-08 20.72670 65.27605 63.59760
## 2007-01-09 20.89846 64.58615 63.53651
tail(data)
## AIR.PA.Adjusted BA.Adjusted LMT.Adjusted
## 2019-03-25 114.04 370.46 294.12
## 2019-03-26 116.34 370.38 294.92
## 2019-03-27 115.48 374.21 293.95
## 2019-03-28 116.26 374.44 295.60
## 2019-03-29 117.90 381.42 300.16
## 2019-04-01 117.76 391.54 304.29
#=================================
# Minimum variance portfolio
#=================================
library(fBasics)
## Loading required package: timeDate
##
## Attaching package: 'timeDate'
## The following objects are masked from 'package:PerformanceAnalytics':
##
## kurtosis, skewness
## Loading required package: timeSeries
##
## Attaching package: 'timeSeries'
## The following object is masked from 'package:zoo':
##
## time<-
##
## Attaching package: 'fBasics'
## The following object is masked from 'package:TTR':
##
## volatility
Sigma = cov(firm_data1[,2:4])
std = sqrt(diag(Sigma))
ones = rep(1,3)
one.vec = matrix(ones, ncol=1)
a = inv(Sigma)%*%one.vec
b = t(one.vec)%*%a
mvp.w =a / as.numeric(b)
mvp.w
## [,1]
## EPS.BA -0.1902673
## EPS.AIR 1.3422878
## EPS.LMT -0.1520205
mvp.ret<-sum((mvp.w)*colMeans(firm_data1[,2:4]))
mvp.ret
## [1] 22.92708
#==================================
# Assume return is specified as 0.06.
# Try to find its optimal weight and standard deviation (tangency portfolio),
# expected return and Sharpe ratio.
#=================================
mu<-0.06/12
return <- firm_data1[,2:4]
Ax <- rbind(2*cov(return), colMeans(return), rep(1, ncol(return)))
Ax <- cbind(Ax, rbind(t(tail(Ax, 2)), matrix(0, 2, 2)))
b0 <- c(rep(0, ncol(return)), mu, 1)
out<-solve(Ax, b0)
wgt<-out[1:3]
wgt
## EPS.BA EPS.AIR EPS.LMT
## -0.2256986 1.4471669 -0.2214683
sum(wgt)
## [1] 1
ret.out<-sum(wgt*colMeans(return))
ret.out.annual<-ret.out*12
ret.out.annual
## [1] 0.06
std.out<-sqrt(t(wgt)%*%cov(return)%*%wgt)
std.out.annual<-std.out*sqrt(12)
std.out.annual
## [,1]
## [1,] 26.38877
#====================================================================================================
# Or write your own function to find min var for a specified return mu;
# Reference: Introduction to R for Quantitative Finance: Chapter 2, p31
# Solve a diverse range of problems with R, one of the most powerful tools for quantitative finance
#====================================================================================================
return = firm_data1[,2:4]
#specified portfolio return: mu
mu=0.06/12
minvariance <- function(return, mu) {
#return <- log(tail(assets, -1) / head(assets, -1))
Ax <- rbind(2*cov(return), colMeans(return), rep(1, ncol(return)))
Ax <- cbind(Ax, rbind(t(tail(Ax, 2)), matrix(0, 2, 2)))
b0 <- c(rep(0, ncol(return)), mu, 1)
zx<-solve(Ax, b0)
weight<-zx[1:ncol(return)]
ret.out<-sum(weight*colMeans(return))
std.out<-sqrt(t(wgt)%*%cov(return)%*%wgt)
list(weight=weight, rtn=ret.out, sd=std.out)
}
minvariance(return, mu)
## $weight
## EPS.BA EPS.AIR EPS.LMT
## -0.2256986 1.4471669 -0.2214683
##
## $rtn
## [1] 0.005
##
## $sd
## [,1]
## [1,] 7.617781
#======================================================
# Create frontier function to plot efficient frontier
#======================================================
frontier <- function(return){
#return <- log(tail(assets, -1) / head(assets, -1))
n = ncol(return)
Q = cov(return)
Ax <- rbind(2*cov(return), colMeans(return), rep(1, n))
Ax <- cbind(Ax, rbind(t(tail(Ax, 2)), matrix(0, 2, 2)))
r <- colMeans(return)
rbase <- seq(min(r), max(r), length = 100)
s <- sapply(rbase, function(x) {
b0 <- c(rep(0, ncol(return)), x, 1)
y <- head(solve(Ax, b0), n)
sqrt(y%*%Q%*%y)
})
plot(s, rbase, xlab = 'Std', ylab = 'Return')
}
frontier(return)
