R Notebook

When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Ctrl+Shift+K to preview the HTML file).

The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike Knit, Preview does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed.S

##Extract Price Series Here we extract the index and the comparison equity We use the Quantmod Libraries standard getSymbols tool. We then pass symbols to Yahoo! Finance via the getSymbols() function from the quantmod package. This will return an object with the opening price, closing price, adjusted price, daily high, daily low and daily volume for each ticker. To isolate the adjusted price, we use the map() function from the purrr package and apply Ad(get(.)) to the imported prices. This will get() the adjusted price from each of our individual price series. If we wanted the closing price, we would run Cl(get(.)). That . refers to our initial object. We could stop here and have the right substance - daily prices for 2 ETFs - but the format would not be great as we would have a list of 2 xts objects. This is because the map() function returns a list by default. The reduce(merge) function will merge the 2 lists into one xts object. The merge() function looks for the date index shared by our objects and uses that index to align the data. After we set prices to daily we estimate returns with “asset_returns_xts <- na.omit(Return.calculate(prices_daily, method =”log“))” na.omit removes blanks, Return.Calculate calculates daily returns, throughout the rest of this we use the dailyReturn() function to achieve the same thing in a simpler fashion, method = log signifies log daily or geometric returns as opposed to the default arithmetic rate of change.

Finally, we want intuitive column names and use colnames<- to rename the columns according to the symbols object. In order to add technical indicators in a chart, Quantmod provides the add{IndicatorName)() function. We will use this function for plotting the Bollinger Band of SPY prices. Both indicators are plotted with the default parameters which are specified below: To add a second series on the same chart it is similar to adding and indicator in this case after you add the first chart series as normal, add this TA = “addTA(i_a.cls, on= 1)”, yrange = c(0, 380) where in the quote addTA in the brackets is the second series, the on = 1 signifies on the same chart not below,and you may need to adjust the y axis as yrange = or; the x axis as xrange =. The performance.Chart shows cumulative daily returns to compare the indices, the in the next band daily range of price returns, then below, the drawdowns, or biggest daily losses.

#startDate<- as.Date("2007-01-03")

ticker <-c("BAC","QQQ","IWM")

# getSymbols(ticker, from = "1993-01-01", auto.assign = TRUE)
# tidyQuant method in timeseries tibble
prices <- 
  getSymbols(ticker, src = 'yahoo', 
             from = "2000-01-01", 
             auto.assign = TRUE, 
             warnings = FALSE) %>% 
  map(~Ad(get(.))) %>% 
  reduce(merge) %>%
  `colnames<-`(ticker)



prices_daily <- to.daily(prices, indexAt = "lastof", OHLC = FALSE)

asset_returns_xts <- na.omit(Return.calculate(prices_daily, method = "log"))


#standard quantmod method
getSymbols("BAC", from= "2000-01-01", auto.assign = TRUE)

getSymbols("IWM", from= "2000-01-01", auto.assign = TRUE)

chartSeries(QQQ)
chartSeries(BAC)
i_a.cls <- BAC$BAC.Adjusted
s_a.cls <- QQQ$QQQ.Adjusted

sht_series <- i_a.cls["2020"]
#sht_series

chartSeries(sht_series)
chartSeries(sht_series,
            theme = chartTheme("white"),
            TA = c(addBBands(n = 20, sd = 2, ma = "SMA", draw = 'bands', on = -1))) 



chartSeries(s_a.cls,TA = "addTA(i_a.cls, on= 1)", yrange = c(0, 380))


charts.PerformanceSummary(asset_returns_xts)

THis small calculation shows how much of the total twenty years of the data the market has been down, it reflects the Hurst Index shown later

n_neg <- asset_returns_xts[asset_returns_xts$BAC<0]


neg <- n_neg$BAC
head(neg)
ndf <- as.data.frame(neg)
head(ndf)

count_neg <- count(ndf)
count_neg
neg_asset <- as.data.frame(asset_returns_xts$BAC)
count_asset <- count(neg_asset)
down <- count_neg/count_asset
down
up <- 1-down
up

Here we calculate a table of Capital Asset Pricing model Metrics such as Sharpe ratio , treynor ratio, beta, and alpha, etc. We also calculate a relative performance chart of IWM against the QQQ

head(asset_returns_xts)




chart.RelativePerformance(asset_returns_xts,asset_returns_xts$QQQ)
table.CAPM(asset_returns_xts,asset_returns_xts$QQQ,Rf = 0, scale = 252)

The probabilities and the expected range are really based on a snapshot in time. Recognition of what your expectations are for implied volatility in the future is therefore a necessary part of assessing the confidence and validity you place in the probability at which you arrive. Many traders do not even look at the historical distribution curve at all, and I think that is a huge mistake. We will talk more about that later, but for now, we will concentrate on our forwardlooking distribution curve as defined by the implied of the stocks. The original stock price model assumed that underlying assets moved randomly. After Merton got involved in the process, it was recognized that the curve was actually lognormal in shape. An asset’s value can theoretically go up infinitely but can never go to a negative value. One last facet of the distribution curve is its level of kurtosis. statistically speaking, a measure of the peakedness of the distribution curve. A statistical distribution where the points along the X‐axis are clustered, resulting in a higher peak (higher kurtosis) than the curvature found in a normal distribution. Th is high peak and corresponding fat tails means the distribution is more clustered around the mean than in a mesokurtic or platykurtic distribution, and will have a relatively smaller standard deviation. As we know, the standard deviation, as described in the definition, is synonymous with volatility. So, when we are in a relatively low volatility environment, the skew takes on a more leptokurtic shape. As a result, the more leptokurtic our distribution curves become, the greater the risk that traders will get caught with their pants down when the underlying has a tail risk–style event.

ts_r <- as.timeSeries(asset_returns_xts$BAC)
tb_r <- as.data.frame(asset_returns_xts$BAC)
histPlot(ts_r)
#chart.Histogram(ts_r,breaks= "FD",methods = c("add.rug","add.risk"),color ="blue", )
chart.Histogram(ts_r,main =" SPY Histogram",methods = c("add.density","add.centered","add.rug","add.risk"),color = bluemono)
library(stargazer)
out <- summary(ts_r)
stargazer(out, type = "text")
table.Stats(ts_r)
#out_rt <- t(table.CalendarReturns(ts_r))
#stargazer(out_rt, type = "text")
#
# layout(rbind(c(1,2),c(3,4)))
#chart.Histogram(asset_returns_xts, main = "Plain", methods = NULL)
 #chart.Histogram(asset_returns_xts, main = "Density", breaks=40,
 #methods = c("add.density", "add.normal"))
 #chart.Histogram(asset_returns_xts, main = "Skew and Kurt", methods = c
 #("add.centered", "add.rug"))
 #chart.Histogram(asset_returns_xts, main = "Risk Measures", methods = c
 #("add.risk"))

The below highlights the relative and cumulative frequency of the data as illustrated in the histograms in a data table. It is usefull for checking, given either the Implied volatility rank or percentile the likelihood a downside breach, using a rule of thumb heuristic based on the table stats.

factorx <- factor(cut(asset_returns_xts$BAC, breaks=nclass.Sturges(asset_returns_xts$BAC)))
#TSPYulate and turn into data.frame
xout <- as.data.frame(table(factorx))
#Add cumFreq and proportions
xout_r <- transform(xout, cumFreq = cumsum(Freq), relative = prop.table(Freq))

xout_r
tb_xout <- tibble(xout_r)


tb_xout

The conditional variance at risk is the variance at risk calculated using actual historical data, whereas the V.A.R. simply uses a z score to estimate risk.

CVaR(asset_returns_xts, alpha = .01, type ="sample")
CVaR(asset_returns_xts, alpha = .05, type ="sample")

Here we load the Pracma and fractal libraries for estimating various forms of the Hurst re scaled range statistic.

library(pracma)
# library(fractal)
library(quantmod)
getSymbols("BAC")

geom_return_BAC <- dailyReturn(BAC, type = "log")

x <- geom_return_BAC

Having estimated daily returns #(geom_return_BAC) and assigned to variable x, we then estimate the Hurst index for the stock at hand. first with Hurstexp(x) from the pracma libraries which gives the simple R/S statistic and some more experimental calculated variations.hurstexp(x) calculates the Hurst exponent of a time series x using R/S analysis, after Hurst, with slightly different approaches, or corrects it with small sample bias, see for example Weron.

These approaches are a corrected R/S method, an empirical and corrected empirical method, and a try at a theoretical Hurst exponent. It should be mentioned that the results are sometimes very different, so providing error estimates will be highly questionable.

Optimal sample sizes are automatically computed with a length that possesses the most divisors among series shorter than x by no more than 1 percent.

Hurstspec() form the fractal library is estimated next, it is estimated by spectral density function. Hurst Coefficient Estimation Via Spectral Regression Function to estimate the Hurst parameter H of a time series by linear regression of the log(spectrum) versus log(frequency) with frequency points accumulated into boxes of equal width on a logarithmic scale and spectrum values averaged over each box.We then look at a variation of this method, then 2 more; a simple user defined R/S H stat function, and a aggabs method of hurstblock function, which calculates greater scaled samples of data for their hurst exponent then averages them.The series is partitioned into m groups. Within each group, the first absolute moment about the mean of the entire series is evaluated. A measure of the variability of this statistic between groups is calculated. The number of groups, m, is increased and the process is repeated. The observed variability changes with increasing m in a way related by theory to the Hurst parameter H of the input series. For the methods used here, a log-log plot of variability versus number of groups is, ideally, linear, with a slope related to H, so H can be determined by linear regression.

hurstexp(x)

# hurstSpec(x)
# 
# hurstSpec(x, method="robinson", sdf.method="multitaper")
# RoverS(geom_return_BAC)
# 
# "aggAbs" <- hurstBlock(x, method="aggAbs", scale.min=8, scale.max=NULL,
#     scale.ratio=8, weight=function(x) rep(1,length(x)), fit=lm)
# 
# plot(aggAbs)


SimpleHurst <- function(y){
      sd.y <- sd(y)
      m <- mean(y)
      y <- y - m
      max.y <- max(cumsum(y))
      min.y <- min(cumsum(y))
      RS <- (max.y - min.y)/sd.y
      H <- log(RS) / log(length(y))
      return(H)
}
SimpleHurst(x)

# Alternate calculation
#hurst(nile.yr.ts)
#hurstexp(nile.yr.ts)

library(fPortfolio)
library(fExtremes)
library(PerformanceAnalytics)
rsk <- table.DownsideRisk(geom_return_BAC,Rf= 0,p = .99, digits = 4, ci = .99)
rsk
CVaR(geom_return_BAC, alpha = .01)
cvarRisk(geom_return_BAC, alpha = .01, weights = 1)
rslt_drw <- table.Drawdowns(geom_return_BAC)
rslt_drw

we will estimate from the law of small numbers the likelihood of left tail or downside events, one the mean likely drop or expected drop, and one the worst case scenario given the data. The next section uses the Peak over threshold method to estimate the extreme value distribution using a pareto curve fitted by maximum likelihood methods. It is often refer to as the law of small numbers, as extreme value theory after the GFC has come to replace VAR as the serious downside risk estimator for quants.

xt <-  -geom_return_BAC

findThreshold(xt)
thrshd <- findThreshold(xt)
thrshd
findThreshold(xt,doplot = TRUE)

fit1 <- fitgpd(xt, thresh = .04, est = "mle")
fit <- fitgpd(xt, thresh = thrshd, est = "mle")








gpdfit <- gpd(xt,threshold =0.01)


names(gpdfit)
par(mfcol = c(2,2))
plot(gpdfit)

we estimate the pareto distribution then we forecast likely daily downside change, using GARCH metods on GARCH(1,1) , and ARMA(1,0) function.

#fit_dta <- geom_return_BAC[geom_return_BAC< .04]





fit_data_garch <- garchFit(formula=~ arma(1, 0) + garch(1, 1),
         data=xt,
         cond.dist="norm", trace=FALSE)

# Produce forecasts
model.forecast <- fGarch::predict(object=fit_data_garch, n.ahead=1)

cat("Forecasts for Extremes")
model.forecast

est <- model.forecast$meanForecast*last(BAC)

est

Here we fit the most likely extreme tail event given the data sample, with two distributions one at the 5% level of confidence and one at the 1% level of cofidence and we then calculate the maximum expected daily drop in an extreme event.

tp <- tailplot(gpdfit)
gpd.q(tp, pp = 0.999, ci.p = 0.95)
quantile(xt, probs = 0.999, type = 1)
gpd.sfall(tp, 0.99)

loss <- .09857*last(BAC$BAC.Adjusted)
loss

E

These are incomplete quantmod getOption functions not currently run.

library(quantmod)
getOptionChain("BAC", Exp = "2021-01")

#your list of expiration dates

Exp=c("2021-01", "2020-01-20")

x=list()
for(i in 1:length(exp)) x<-c(getOptionChain("BAC",exp[i]),x)

x

E