Load Libraries

library(quantmod)
library(tidyquant)
library(forecast)
library(fBasics)
library(PerformanceAnalytics)
library(tidyverse)
library(tsbox)
library(moments)

Get SPY and SPDR XLE Share prices

c for combining symbols using QUANTMOD package getsymbols command to down load data, arrange using TIDYVERSE package dplyr, date , stock symbol, and adjusted prices.

Turn them into monthly data, with PERFORMANCEANALYTICS package library command to.monthly.

Then head command to check two lines of data.

symbols <- c("SPY","TPR")

prices <- getSymbols(symbols,
             src = 'yahoo',
             from = "2010-01-01",
             auto.assign = TRUE,
             warnings = FALSE) %>%
  map(~Ad(get(.))) %>%
  reduce(merge) %>%
  `colnames<-`(symbols)



prices_monthly <- to.monthly(prices,
                             indexAt = "lastof",
                             OHLC = FALSE)
head(prices_monthly, 3)
##                 SPY      TPR
## 2010-01-31 84.91402 25.26676
## 2010-02-28 87.56292 26.39681
## 2010-03-31 92.89368 28.68692

Calculate daily price changes as log returns.

Use QUANTMOD dailyReturn command to estimate daily returns,log = logarithmic returns (natural logarithm).

Then calculate monthly returns with performanceAnalytics command Return.Calculate.

Check data with head command.

 p <- dailyReturn(TPR, type = "log")
#  fr <- dailyReturn(F, type = "log")
# bacr <- dailyReturn(BAC, type = "log")
# tr <- dailyReturn(T, type = "log")
stk_rtns <-
  Return.calculate(prices_monthly,
                   method = "log") %>%
  na.omit()
head(stk_rtns, 3)
##                   SPY        TPR
## 2010-02-28 0.03071844 0.04375328
## 2010-03-31 0.05909796 0.08319836
## 2010-04-30 0.01535180 0.05489292
head(p)
##            daily.returns
## 2010-01-04 -0.0158472715
## 2010-01-05  0.0120450127
## 2010-01-06  0.0194024138
## 2010-01-07  0.0005336446
## 2010-01-08 -0.0058855706
## 2010-01-11  0.0277832396

Show Charts, histogram and Conditional Variance at Risk.

First for Rmetrics histPlot the xts data must be converted to timeSeries format.

Then CVAR is or Conditional variance at risk which is the historical datas variance at risk at either 99% or 95% in this instance 95%.

Then we fit an autoArima model from the FORECAST package, just out of curiousity to see whether any autocorrelation or trends exist. If so I t would take in depth analysis beyond this brief intro .

#timeSeries conversion for Rmetrics libraries

p_ts<- as.timeSeries(p)
#ts_T <- as.timeSeries(tr)
#tS_BAC <- as.timeSeries(bacr)
#fbasics Histogram...histPlot
histPlot(p_ts,col ="steelblue",fit=TRUE,title=TRUE,grid=TRUE,rug=TRUE)

# histPlot(tS_T,col ="steelblue",fit=TRUE,title=TRUE,grid=TRUE,rug=TRUE)
# histPlot(tS_BAC,col ="steelblue",fit=TRUE,title=TRUE,grid=TRUE,rug=TRUE)
#peXLEormanceAnalytics histogram 
chart.Histogram(p_ts,methods=c("add.risk","add.qqplot"))

#peXLEormanceAnalytics Conditional Variance at Risk "CVAR"
CVaR(p_ts)
##    daily.returns
## ES    -0.0832222
head(p_ts)
## GMT
##            daily.returns
## 2010-01-04 -0.0158472715
## 2010-01-05  0.0120450127
## 2010-01-06  0.0194024138
## 2010-01-07  0.0005336446
## 2010-01-08 -0.0058855706
## 2010-01-11  0.0277832396
#select range for Arima plot/forecat
# Z <- spr["2020/2022"]
# Z <- as.ts(Z)
# fit <- auto.arima(Z)
# autoplot(Z)
# ts.plot(Z)
# fit
# ffit <- forecast(fit, h = 12)
# autoplot(ffit)
# Acf(Z)

Performance Analytics analysis of data under question.

We look at a comparative cumulative performance chart PerformanceSummary.

Then annual and calendar returns.

Then a table of Capital asset Portfolio management benchmark statistics.

#cumulative returns, monthly returns & drawdowns in peXLEormanceAnalytics summary chart.
charts.PerformanceSummary(stk_rtns,TPR= 0,SPY)

#peXLEormanceAnalytics tables.
table.AnnualizedReturns(stk_rtns$TPR)
table.CalendarReturns(stk_rtns$TPR)
table.CAPM(stk_rtns$TPR,stk_rtns$SPY)
table.AnnualizedReturns(stk_rtns$SPY)
table.CalendarReturns(stk_rtns$SPY)
#table.AnnualizedReturns(stk_rtns$QQQ)
#table.CalendarReturns(stk_rtns$QQQ)

Here we compare skewness and Kurtosis.

Skewness is the tilt or spread to either positive or negative tendency among the data.

Kurtosis is the level at which the peak or central tendency of the data reaches and the width of it spread or the extent to which it flattens towards ’ fat tails’ or not.

# kurtosis(xlr)
# skewness(xlr)
# SkewnessKurtosisRatio(xlr)
kurtosis(TPR$TPR.Adjusted)
## TPR.Adjusted 
##     3.080122
skewness(TPR$TPR.Adjusted)
## TPR.Adjusted 
##   -0.2713978
SkewnessKurtosisRatio(TPR$TPR.Adjusted)
## [1] -0.08811268
kurtosis(SPY$SPY.Adjusted)
## SPY.Adjusted 
##     2.916123
skewness(SPY$SPY.Adjusted)
## SPY.Adjusted 
##    0.8386827
SkewnessKurtosisRatio(SPY$SPY.Adjusted)
## [1] 0.287602

Here we have a frequency table in a tibble format.

factorx <- factor(cut(stk_rtns$TPR, breaks=nclass.Sturges(stk_rtns$TPR)))
#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
# factorx <- factor(cut(stk_rtns$QQQ, breaks=nclass.Sturges(stk_rtns$XLE)))
# #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

Here we have the Hurst exponent

Here we estimate the hurst exponent with a manual function and a variety of estimates from the pracma package library. It describes the long run tendency of the data to mean revert or trend( either direction).

#p_ts<- as.timeSeries(p)

#install.packages(pracma)
library(pracma)

hurstexp(p)
## Simple R/S Hurst estimation:         0.5448802 
## Corrected R over S Hurst exponent:   0.5670491 
## Empirical Hurst exponent:            0.5534556 
## Corrected empirical Hurst exponent:  0.5183973 
## Theoretical Hurst exponent:          0.532951
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(p)
## [1] 0.5360987

Here we estimate the likelihood of a changepoint with bayesian changepoint package library functions.

library(bcp)
## Warning: package 'bcp' was built under R version 4.0.5
## Loading required package: grid
y <- stk_rtns$TPR
out <- bcp(y, burnin=500, mcmc=500)
str(y)
## An 'xts' object on 2010-02-28/2022-05-31 containing:
##   Data: num [1:148, 1] 0.0438 0.0832 0.0549 -0.0154 -0.114 ...
##  - attr(*, "dimnames")=List of 2
##   ..$ : NULL
##   ..$ : chr "TPR"
##   Indexed by objects of class: [Date] TZ: UTC
##   xts Attributes:  
## List of 5
##  $ src             : chr "yahoo"
##  $ updated         : POSIXct[1:1], format: "2022-05-06 07:03:05"
##  $ ret_type        : chr "log"
##  $ coredata_content: chr "logReturn"
##  $ na.action       : 'omit' int 1
##   ..- attr(*, "index")= num 1.26e+09
#library(bcp)
fit_bcp = bcp(y, d = 10)
plot(fit_bcp)

x <- stk_rtns$SPY
out <- bcp(x, burnin=500, mcmc=500)
str(x)
## An 'xts' object on 2010-02-28/2022-05-31 containing:
##   Data: num [1:148, 1] 0.0307 0.0591 0.0154 -0.0828 -0.0531 ...
##  - attr(*, "dimnames")=List of 2
##   ..$ : NULL
##   ..$ : chr "SPY"
##   Indexed by objects of class: [Date] TZ: UTC
##   xts Attributes:  
## List of 5
##  $ src             : chr "yahoo"
##  $ updated         : POSIXct[1:1], format: "2022-05-06 07:03:05"
##  $ ret_type        : chr "log"
##  $ coredata_content: chr "logReturn"
##  $ na.action       : 'omit' int 1
##   ..- attr(*, "index")= num 1.26e+09
#library(bcp)
fit_bcp = bcp(x, d = 10)
plot(fit_bcp)

chartSeries(TPR$TPR.Adjusted)

Here we estimate the extreme distribution of tail data with a general Pareto distribution the Peak over Threshold extreme risk. This is an absolute robust risk measure, as opposed to VAR and CVAR and monte carlo VAR

library(quantmod)
library(evir)
library(fExtremes)


#getSymbols("TPR", from = "2010-01-01")

# s <- dailyReturn(TPR, type = "log")

s <- p

emplot(s)

meplot(s)

fit <- gpd(s, threshold = NA, nextremes = 20)
fit
## $n
## [1] 3107
## 
## $data
##  [1] 0.11257363 0.08271187 0.09352957 0.09366143 0.10780975 0.08623276
##  [7] 0.11341900 0.08146761 0.13885258 0.25988920 0.08250659 0.16197776
## [13] 0.10405598 0.09837134 0.09307221 0.07887492 0.08035921 0.08142304
## [19] 0.08045781 0.10109614
## 
## $threshold
## [1] 0.07814474
## 
## $p.less.thresh
## [1] 0.9935629
## 
## $n.exceed
## [1] 20
## 
## $method
## [1] "ml"
## 
## $par.ests
##         xi       beta 
## 0.41004923 0.01751092 
## 
## $par.ses
##         xi       beta 
## 0.32934781 0.00673566 
## 
## $varcov
##              [,1]          [,2]
## [1,]  0.108469980 -1.382439e-03
## [2,] -0.001382439  4.536912e-05
## 
## $information
## [1] "observed"
## 
## $converged
## [1] 0
## 
## $nllh.final
## [1] -52.69487
## 
## attr(,"class")
## [1] "gpd"
# run the plot below it has an interactive menu so I can not knit the document it freezes #waiting for menu response.
#plot(fit)

# 1
# 2
# 3
# 4
# 5
# 0

tp <-  tailplot(fit)



gpd.q(tp ,pp = .999, ci.p = .95)
##  Lower CI  Estimate  Upper CI 
## 0.1056377 0.1270782 0.2023230
gpd.sfall(tp, .99)

##   Lower CI   Estimate   Upper CI 
## 0.08775067 0.09586440 0.38983379

Here we estimate Bollinger bands.

chartSeries(SPY$SPY.Adjusted, TA="addBBands(n=20)", subset="2020-12-01::2022-03-22")

chartSeries(TPR$TPR.Adjusted, TA="addBBands(n=20)", subset="2020-12-01::2022-03-22")