Load Libraries

Get 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("RIO.AX","BHP.AX","OZL.AX","NCM.AX","NST.AX")

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





prices_monthly <- to.monthly(prices,
                             indexAt = "lastof",
                             OHLC = FALSE)
head(prices_monthly, 3)
##              RIO.AX   BHP.AX   OZL.AX   NCM.AX    NST.AX
## 2019-04-30 85.72207 30.19167 9.462736 23.75175  7.788488
## 2019-05-31 90.15333 30.48228 8.569667 25.84107  9.233970
## 2019-06-30 93.26331 33.22698 9.529242 30.34260 11.078862
stk_rtns_m <-
  Return.calculate(prices_monthly,
                   method = "log") %>%
  na.omit()

head(stk_rtns_m)
##                 RIO.AX       BHP.AX      OZL.AX      NCM.AX      NST.AX
## 2019-05-31  0.05040156  0.009579637 -0.09913268  0.08430883  0.17024233
## 2019-06-30  0.03391487  0.086216507  0.10613630  0.16058757  0.18214989
## 2019-07-31 -0.04797139 -0.009765717  0.02559192  0.10761152  0.11041224
## 2019-08-31 -0.11161906 -0.116158914 -0.11741875  0.04528160 -0.08918001
## 2019-09-30  0.05649216  0.043585510  0.05783195 -0.06273831 -0.07501329
## 2019-10-31 -0.01950491 -0.021192434  0.05466958 -0.10456141 -0.11244950

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(prices$RIO.AX, type = "log")
   q <- dailyReturn(prices$OZL.AX, type = "log")
#  fr <- dailyReturn(F, type = "log")
# bacr <- dailyReturn(BAC, type = "log")
# tr <- dailyReturn(T, type = "log")

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 .

##    daily.returns
## ES   -0.04541402
## Series: p 
## ARIMA(0,0,0) with non-zero mean 
## 
## Coefficients:
##        mean
##       1e-04
## s.e.  7e-04
## 
## sigma^2 estimated as 0.000355:  log likelihood=2116.72
## AIC=-4229.44   AICc=-4229.43   BIC=-4220

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_m, RF= 0)

#peXLEormanceAnalytics tables.
table.AnnualizedReturns(stk_rtns_m$RIO.AX)
##                           RIO.AX
## Annualized Return         0.0001
## Annualized Std Dev        0.2602
## Annualized Sharpe (Rf=0%) 0.0005
table.CalendarReturns(stk_rtns_m$OZL.AX)
##        Jan   Feb   Mar  Apr  May   Jun  Jul   Aug  Sep  Oct Nov  Dec OZL.AX
## 2019    NA    NA    NA   NA -9.9  10.6  2.6 -11.7  5.8  5.5 3.9  0.0    4.5
## 2020  -3.8 -12.1 -18.4 19.9  4.9  15.1 21.8   7.0 -3.7  5.9 9.4 14.6   66.3
## 2021  -1.2  18.3   2.8  4.7  5.5 -11.8  3.1   2.1 -4.6 11.0 3.4  8.2   45.8
## 2022 -15.5   6.0   4.5 -6.5 -3.7 -30.5 -2.5    NA   NA   NA  NA   NA  -42.8
#table.CAPM(stk_rtns_m$SFY,stk_rtns_m$RIO.AX)
table.AnnualizedReturns(stk_rtns_m$RIO.AX)
##                           RIO.AX
## Annualized Return         0.0001
## Annualized Std Dev        0.2602
## Annualized Sharpe (Rf=0%) 0.0005
table.CalendarReturns(stk_rtns_m$RIO.AX)
##       Jan   Feb   Mar  Apr May   Jun  Jul   Aug   Sep   Oct Nov  Dec RIO.AX
## 2019   NA    NA    NA   NA 5.0   3.4 -4.8 -11.2   5.6  -2.0 6.4  3.5    4.8
## 2020 -1.6 -12.4   0.8  3.4 6.5   4.8  4.0  -1.9  -3.8  -2.0 9.3 11.6   17.6
## 2021 -3.1  14.2 -12.9  9.0 2.1   2.3  5.2 -15.5 -11.2 -10.4 3.5  6.8  -14.2
## 2022 10.8   5.8   1.5 -5.4 1.4 -10.8 -7.1    NA    NA    NA  NA   NA   -5.4
#table.AnnualizedReturns(stk_rtns$QQQ)
#table.CalendarReturns(stk_rtns$QQQ)
table.Stats(stk_rtns_m$RIO.AX)
##                  RIO.AX
## Observations    39.0000
## NAs              0.0000
## Minimum         -0.1549
## Quartile 1      -0.0431
## Median           0.0211
## Arithmetic Mean  0.0028
## Geometric Mean   0.0000
## Quartile 3       0.0543
## Maximum          0.1424
## SE Mean          0.0120
## LCL Mean (0.95) -0.0215
## UCL Mean (0.95)  0.0272
## Variance         0.0056
## Stdev            0.0751
## Skewness        -0.4040
## Kurtosis        -0.6804

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(p)
## daily.returns 
##      4.483183
skewness(p)
## daily.returns 
##    -0.2099957
SkewnessKurtosisRatio(p)
## [1] -0.04684076
kurtosis(q)
## daily.returns 
##      6.205117
skewness(q)
## daily.returns 
##    -0.5262998
SkewnessKurtosisRatio(q)
## [1] -0.08481707

Here we have a frequency table in a tibble format.

factorx <- factor(cut(p, breaks=nclass.Sturges(p)))
#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
##               factorx Freq cumFreq    relative
## 1    (-0.0765,-0.062]    4       4 0.004825090
## 2    (-0.062,-0.0475]    8      12 0.009650181
## 3   (-0.0475,-0.0331]   25      37 0.030156815
## 4   (-0.0331,-0.0187]   76     113 0.091676719
## 5  (-0.0187,-0.00432]  207     320 0.249698432
## 6   (-0.00432,0.0101]  282     602 0.340168878
## 7     (0.0101,0.0245]  169     771 0.203860072
## 8     (0.0245,0.0389]   42     813 0.050663450
## 9     (0.0389,0.0533]   12     825 0.014475271
## 10    (0.0533,0.0677]    3     828 0.003618818
## 11    (0.0677,0.0823]    1     829 0.001206273
tb_xout <- tibble(xout_r)


tb_xout
## # A tibble: 11 x 4
##    factorx             Freq cumFreq relative
##    <fct>              <int>   <int>    <dbl>
##  1 (-0.0765,-0.062]       4       4  0.00483
##  2 (-0.062,-0.0475]       8      12  0.00965
##  3 (-0.0475,-0.0331]     25      37  0.0302 
##  4 (-0.0331,-0.0187]     76     113  0.0917 
##  5 (-0.0187,-0.00432]   207     320  0.250  
##  6 (-0.00432,0.0101]    282     602  0.340  
##  7 (0.0101,0.0245]      169     771  0.204  
##  8 (0.0245,0.0389]       42     813  0.0507 
##  9 (0.0389,0.0533]       12     825  0.0145 
## 10 (0.0533,0.0677]        3     828  0.00362
## 11 (0.0677,0.0823]        1     829  0.00121
factorx <- factor(cut(q, breaks=nclass.Sturges(q)))
#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
##                factorx Freq cumFreq    relative
## 1      (-0.153,-0.131]    1       1 0.001206273
## 2      (-0.131,-0.109]    1       2 0.001206273
## 3     (-0.109,-0.0875]    2       4 0.002412545
## 4    (-0.0875,-0.0658]    2       6 0.002412545
## 5    (-0.0658,-0.0441]   25      31 0.030156815
## 6    (-0.0441,-0.0224]   79     110 0.095295537
## 7  (-0.0224,-0.000659]  276     386 0.332931242
## 8   (-0.000659,0.0211]  299     685 0.360675513
## 9      (0.0211,0.0428]  117     802 0.141133896
## 10     (0.0428,0.0645]   21     823 0.025331725
## 11     (0.0645,0.0864]    6     829 0.007237636
tb_xout <- tibble(xout_r)


tb_xout
## # A tibble: 11 x 4
##    factorx              Freq cumFreq relative
##    <fct>               <int>   <int>    <dbl>
##  1 (-0.153,-0.131]         1       1  0.00121
##  2 (-0.131,-0.109]         1       2  0.00121
##  3 (-0.109,-0.0875]        2       4  0.00241
##  4 (-0.0875,-0.0658]       2       6  0.00241
##  5 (-0.0658,-0.0441]      25      31  0.0302 
##  6 (-0.0441,-0.0224]      79     110  0.0953 
##  7 (-0.0224,-0.000659]   276     386  0.333  
##  8 (-0.000659,0.0211]    299     685  0.361  
##  9 (0.0211,0.0428]       117     802  0.141  
## 10 (0.0428,0.0645]        21     823  0.0253 
## 11 (0.0645,0.0864]         6     829  0.00724

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.5016713 
## Corrected R over S Hurst exponent:   0.5018142 
## Empirical Hurst exponent:            0.4140173 
## Corrected empirical Hurst exponent:  0.3828936 
## Theoretical Hurst exponent:          0.5377428
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.504073

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 <- p
out <- bcp(y, burnin=500, mcmc=500)
str(y)
## An 'xts' object on 2019-04-12/2022-07-21 containing:
##   Data: num [1:829, 1] 0 0.00248 0.00446 -0.04859 0.0138 ...
##  - attr(*, "dimnames")=List of 2
##   ..$ : NULL
##   ..$ : chr "daily.returns"
##   Indexed by objects of class: [Date] TZ: UTC
##   xts Attributes:  
##  NULL
#library(bcp)
fit_bcp = bcp(y, d = 10)
plot(fit_bcp)

x <- q
out <- bcp(x, burnin=500, mcmc=500)
str(x)
## An 'xts' object on 2019-04-12/2022-07-21 containing:
##   Data: num [1:829, 1] 0 0.00569 -0.00284 -0.00761 0.00571 ...
##  - attr(*, "dimnames")=List of 2
##   ..$ : NULL
##   ..$ : chr "daily.returns"
##   Indexed by objects of class: [Date] TZ: UTC
##   xts Attributes:  
##  NULL
#library(bcp)
fit_bcp = bcp(x, d = 10)
plot(fit_bcp)

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("GOLD", from = "2010-01-01")

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

s <- p

emplot(s)

meplot(s)

fit <- gpd(s, threshold = NA, nextremes = 20)
fit
## $n
## [1] 829
## 
## $data
##  [1] 0.04644962 0.06662906 0.05228578 0.05014134 0.04429245 0.05625985
##  [7] 0.04567389 0.04018803 0.03617713 0.03802228 0.06641780 0.08214717
## [13] 0.03570095 0.04484499 0.03920479 0.04046767 0.04008610 0.04513794
## [19] 0.03656172 0.03999480
## 
## $threshold
## [1] 0.03545233
## 
## $p.less.thresh
## [1] 0.9758745
## 
## $n.exceed
## [1] 20
## 
## $method
## [1] "ml"
## 
## $par.ests
##           xi         beta 
## -0.001242932  0.011897114 
## 
## $par.ses
##          xi        beta 
## 0.273702319 0.004113958 
## 
## $varcov
##               [,1]          [,2]
## [1,]  0.0749129596 -8.724759e-04
## [2,] -0.0008724759  1.692465e-05
## 
## $information
## [1] "observed"
## 
## $converged
## [1] 0
## 
## $nllh.final
## [1] -68.6549
## 
## 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.06075377 0.07324920 0.12322076
gpd.sfall(tp, .99)

##   Lower CI   Estimate   Upper CI 
## 0.05034444 0.05779352 0.11785092

Here we estimate Bollinger bands.

chartSeries(p$daily.returns,theme= "white", TA="addBBands(n=20)", subset="2020-12-01::2022-03-22")

chartSeries(prices$RIO.AX,theme = "white", TA="addBBands(n=20)", subset="2020-12-01::2022-03-22")