library(quantmod)
library(tidyquant)
library(forecast)
library(fBasics)
library(PerformanceAnalytics)
library(tidyverse)
symbols <- c("SPY","HPE")
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)
## Warning in to.period(x, "months", indexAt = indexAt, name = name, ...): missing
## values removed from data
head(prices_monthly, 3)
## SPY HPE
## 2015-10-31 185.8525 7.264253
## 2015-11-30 186.5318 7.333339
## 2015-12-31 183.3081 7.528334
spr <- dailyReturn(prices$SPY, type = "log")
hpr <- dailyReturn(prices$HPE, type = "log")
stk_rtns <-
Return.calculate(prices_monthly,
method = "log") %>%
na.omit()
head(stk_rtns, 3)
## SPY HPE
## 2015-11-30 0.003648498 0.009465467
## 2015-12-31 -0.017433686 0.026242832
## 2016-01-31 -0.051068570 -0.099529802
#timeSeries conversion for Rmetrics libraries
tS_HPE <- as.timeSeries(hpr)
#fbasics Histogram...histPlot
histPlot(tS_HPE,col ="steelblue",fit=TRUE,title=TRUE,grid=TRUE,rug=TRUE)
#performanceAnalytics histogram
chart.Histogram(hpr,methods=c("add.risk","add.qqplot"))
#performanceAnalytics Conditional Variance at Risk "CVAR"
CVaR(tS_HPE)
## daily.returns
## ES -0.07440194
head(prices$HPE)
## HPE
## 2010-01-04 NA
## 2010-01-05 NA
## 2010-01-06 NA
## 2010-01-07 NA
## 2010-01-08 NA
## 2010-01-11 NA
#select range for Arima plot/forecat
hpe <- prices$HPE["2016/2021"]
hp <- as.ts(hpe)
fit <- auto.arima(hp)
autoplot(hp)
fit
## Series: hp
## ARIMA(0,1,0)
##
## sigma^2 estimated as 0.0639: log likelihood=-65.65
## AIC=133.29 AICc=133.29 BIC=138.61
#cumulative returns, monthly returns & drawdowns in performanceAnalytics summary chart.
charts.PerformanceSummary(stk_rtns,Rf = 0)
#performanceAnalytics tables.
table.AnnualizedReturns(stk_rtns$HPE)
table.CalendarReturns(stk_rtns$HPE)
table.CAPM(stk_rtns$HPE,stk_rtns$SPY)
table.AnnualizedReturns(stk_rtns$SPY)
table.CalendarReturns(stk_rtns$SPY)
factorx <- factor(cut(stk_rtns$HPE, breaks=nclass.Sturges(stk_rtns$HPE)))
#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
ts_hp <- as.timeSeries(stk_rtns$HPE)
library(pracma)
hurstexp(ts_hp)
## Simple R/S Hurst estimation: 0.5655559
## Corrected R over S Hurst exponent: 0.6739064
## Empirical Hurst exponent: 0.5230037
## Corrected empirical Hurst exponent: 0.5664201
## Theoretical Hurst exponent: 0.5184482
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(ts_hp)
## [1] 0.5655559
library(bcp)
## Warning: package 'bcp' was built under R version 4.0.5
## Loading required package: grid
y <- stk_rtns$HPE
out <- bcp(y, burnin=500, mcmc=500)
str(y)
## An 'xts' object on 2015-11-30/2021-12-31 containing:
## Data: num [1:74, 1] 0.00947 0.02624 -0.09953 -0.03626 0.29332 ...
## - attr(*, "dimnames")=List of 2
## ..$ : NULL
## ..$ : chr "HPE"
## Indexed by objects of class: [Date] TZ: UTC
## xts Attributes:
## List of 5
## $ src : chr "yahoo"
## $ updated : POSIXct[1:1], format: "2021-12-18 10:28:51"
## $ na.action : 'omit' int 1
## ..- attr(*, "index")= num 1.45e+09
## $ ret_type : chr "log"
## $ coredata_content: chr "logReturn"
#library(bcp)
fit_bcp = bcp(y, d = 10)
plot(fit_bcp)
chartSeries(prices$HPE)
library(quantmod)
library(evir)
library(fExtremes)
getSymbols("HPE")
## [1] "HPE"
s <- dailyReturn(HPE, type = "log")
emplot(s)
meplot(s)
fit <- gpd(s, threshold = NA, nextremes = 20)
fit
## $n
## [1] 1554
##
## $data
## [1] 0.08096907 0.05869824 0.07301370 0.12689175 0.06549962 0.10022848
## [7] 0.06245995 0.07801636 0.06578335 0.13257308 0.11458303 0.05542453
## [13] 0.05733093 0.06760011 0.09221188 0.05837692 0.05709796 0.05572472
## [19] 0.06097922 0.06724557
##
## $threshold
## [1] 0.05445134
##
## $p.less.thresh
## [1] 0.98713
##
## $n.exceed
## [1] 20
##
## $method
## [1] "ml"
##
## $par.ests
## xi beta
## 0.15161075 0.01889054
##
## $par.ses
## xi beta
## 0.357797795 0.007886003
##
## $varcov
## [,1] [,2]
## [1,] 0.128019262 -2.235566e-03
## [2,] -0.002235566 6.218904e-05
##
## $information
## [1] "observed"
##
## $converged
## [1] 0
##
## $nllh.final
## [1] -56.35677
##
## 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.09102221 0.11339568 0.19885962
gpd.sfall(tp, .99)
## Lower CI Estimate Upper CI
## 0.07190327 0.08244468 0.19885962