Packages

pacman::p_load(dplyr, tidyr, EnvStats, ggplot2, comprehenr, fGarch, LaplacesDemon,
               plotly,laeken, readxl, readr, psych, lattice, GGally, ggpubr, tseries, 
               urca, TSA,lubridate, evir, eva,
               data.table, boot, pharmr, coin, QRM, stats, copula, spatstat, ismev, torch)
year.ret <- function(ts, period = 1) {
  data = c()
  count = 1
  for (i in seq(1, length(ts) - period*52)) {
    if (period == 1) {
      data[count] <- (1 + ts[i:(i+period*52-1)]) %>% prod() - 1
      count = count + 1
    }
    else {
      data[count] <- ((1 + ts[i:(i+period*52-1)]) %>% prod())^(1/period) - 1
      count = count + 1
    }
  }
  return(data)
}
all_data <-  read_csv("DynAll/all_data2.csv")[,c(4, 5, 2, 3, 6)]
## Rows: 1281 Columns: 7
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## dbl  (6): RusEquity_return, RusBonds_return, Gold_return, Gold_USD_return, U...
## date (1): Date
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
#all_data <- all_data[, 2:length(all_data)]
all_data <- all_data %>% replace(all_data == 0, NaN)

#add_data <- read_csv("DynAll/divd_moex_bond.csv")
all_data_cop <- read_csv("DynAll/all_data_monthly2.csv")[,c(4, 5, 2, 3, 6)]
## Rows: 293 Columns: 6
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## dbl  (5): RusEquity_return, RusBonds_return, Gold_USD_return, USD_return, in...
## date (1): Date
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
copt2 <- fit.tcopula(all_data_cop %>% drop_na(), method = "all")
t.cop <- tCopula(
  P2p(copt2$P), 
  dim = 5, 
  dispstr = "un", 
  df = 4)

x <- rCopula(100000, t.cop)

copt2$P
##             [,1]        [,2]        [,3]        [,4]        [,5]
## [1,]  1.00000000 -0.17922842  0.13534448  0.05643076  0.10713495
## [2,] -0.17922842  1.00000000 -0.09828003 -0.21217940 -0.15822266
## [3,]  0.13534448 -0.09828003  1.00000000  0.27116477 -0.03411631
## [4,]  0.05643076 -0.21217940  0.27116477  1.00000000 -0.07002920
## [5,]  0.10713495 -0.15822266 -0.03411631 -0.07002920  1.00000000
all_data %>% drop_na() %>% cor()
##                  Gold_return Gold_USD_return RusEquity_return RusBonds_return
## Gold_return       1.00000000      0.67233415       0.04379932     -0.05742696
## Gold_USD_return   0.67233415      1.00000000       0.19762326      0.09615949
## RusEquity_return  0.04379932      0.19762326       1.00000000      0.26418438
## RusBonds_return  -0.05742696      0.09615949       0.26418438      1.00000000
## USD_return        0.64715627     -0.12886925      -0.14504625     -0.18072901
##                  USD_return
## Gold_return       0.6471563
## Gold_USD_return  -0.1288692
## RusEquity_return -0.1450462
## RusBonds_return  -0.1807290
## USD_return        1.0000000
x %>% cor()
##             [,1]        [,2]        [,3]        [,4]        [,5]
## [1,]  1.00000000 -0.16556989  0.12064037  0.05042423  0.09866808
## [2,] -0.16556989  1.00000000 -0.09038438 -0.19648459 -0.14238535
## [3,]  0.12064037 -0.09038438  1.00000000  0.25268633 -0.02873913
## [4,]  0.05042423 -0.19648459  0.25268633  1.00000000 -0.07326493
## [5,]  0.09866808 -0.14238535 -0.02873913 -0.07326493  1.00000000

##——-АКЦИИ——–

ass_num = 3
asset = unlist(all_data[, ass_num] %>% drop_na()) %>% as.numeric
asset.log = log(asset - min(asset))
dens = density(asset.log, kernel = "optcosine", adjust = 4)
cdff = CDF(dens)
invcdf <- function(q) {
  uniroot(function(x) {cdff(x) - q}, range(asset.log), lower = -100)$root
}
icdf.equity = to_vec(for (i in x[,ass_num]) exp(as.numeric(try(invcdf(i), silent = TRUE))) + min(asset))


#dens = density(asset, kernel = "optcosine", adjust = 3.5)
#cdff = CDF(dens)
#invcdf <- function(q) {
#  uniroot(function(x) {cdff(x) - q}, range(asset), lower = -100)$root
#}
#icdf.equity = to_vec(for (i in x[,ass_num]) as.numeric(try(invcdf(i), silent = TRUE)))
p = 0.05

pareto.lower = gpdFit(-asset, threshold = p, method = "mle")
pareto.upper = gpdFit(asset, threshold = p, method = "mle")

upper.loc = exp(as.numeric(try(invcdf(1-p), silent = TRUE))) + min(asset)
lower.loc = -(exp(as.numeric(try(invcdf(p), silent = TRUE))) + min(asset))

#upper.loc = as.numeric(try(invcdf(1-p), silent = TRUE))
#lower.loc = -as.numeric(try(invcdf(p), silent = TRUE))

equity = to_vec(
  for (i in x[,ass_num]) 
    if (i >= (1 - p)) as.numeric(try(qgpd((i - (1 - p))/(max(x[,ass_num][x[,ass_num] > (1 - p)]) 
                                                         - min(x[,ass_num][x[,ass_num] > (1 - p)])), 
                                          loc = upper.loc, 
                                          scale = pareto.upper$par.ests[1], 
                                          shape = pareto.upper$par.ests[2]), silent = TRUE))
    else if (i <= p) as.numeric(try(-qgpd(i/(max(x[,ass_num][x[,ass_num] <= p]) - 
                                               min(x[,ass_num][x[,ass_num] <= p])), 
                                            loc = lower.loc, 
                                            scale = pareto.lower$par.ests[1], 
                                            shape = pareto.lower$par.ests[2]), silent = TRUE))
    else exp(as.numeric(try(invcdf(i), silent = TRUE))) + min(asset)
    #else as.numeric(try(invcdf(i), silent = TRUE))
  )

hist(asset, breaks = 100, col = rgb(0, 0, 1, 1/4), freq = F)
hist(icdf.equity, breaks = 100, add = T, col = rgb(1, 0, 0, 1/4), freq = F)
hist(equity,
     breaks = 200,
     add = T,
     col = rgb(0, 1, 0, 1/4),
     freq = F)

asset %>% median(na.rm = T)
## [1] 0.002767455
icdf.equity %>% median(na.rm = T)
## [1] 0.002218076
equity %>% median(na.rm = T)
## [1] 0.002227519
y1 <- year.ret(equity + 0.0007, 1)
y1 %>% median(na.rm = T)
## [1] 0.1835578
y1 %>% quantile(na.rm = T, probs = c(0.001, 0.01, .1, 0.9, 0.99, .999))
##       0.1%         1%        10%        90%        99%      99.9% 
## -0.4104333 -0.3107949 -0.1195031  0.5631974  0.9713716  1.3147052
y1.log <- year.ret(icdf.equity + 0.0007, 1)
y1.log %>% median(na.rm = T)
## [1] 0.1895159
y1.log %>% quantile(na.rm = T, probs = c(0.001, 0.01, .1, 0.9, 0.99, .999))
##       0.1%         1%        10%        90%        99%      99.9% 
## -0.4086844 -0.2963056 -0.1109625  0.5782416  0.9929960  1.4218533
y1.h <- year.ret(asset, 1)
y1.h %>% median(na.rm = T)
## [1] 0.1183681
y1.h %>% quantile(na.rm = T, probs = c(0.001, 0.01, .1, 0.9, 0.99, .999))
##        0.1%          1%         10%         90%         99%       99.9% 
## -0.50969170 -0.45427811 -0.08639247  0.48860993  1.25184464  1.51728711
hist(y1.h, 
     freq = F, breaks = 100, col = rgb(1, 0, 0, 1/4))
hist(y1, freq = F, breaks = 100, add = T, col = rgb(0, 1, 0, 1/4))

hist(y1.h, 
     freq = F, breaks = 100, col = rgb(1, 0, 0, 1/4))
hist(y1.log, freq = F, breaks = 100, add = T, col = rgb(0, 1, 0, 1/4))

ecdfPlot(
  equity[equity > quantile(equity, probs = .95, na.rm = T)] - 
         quantile(equity, probs = .95, na.rm = T), ecdf.col = rgb(1, 0, 0, 1/4))
ecdfPlot(
  asset[asset > quantile(asset, probs = .95, na.rm = T)] - 
         quantile(asset, probs = .95, na.rm = T), add = T, ecdf.col = rgb(0, 1, 0, 1/4))

ecdfPlot(
  icdf.equity[icdf.equity > quantile(icdf.equity, probs = .95, na.rm = T)] - 
         quantile(icdf.equity, probs = .95, na.rm = T), ecdf.col = rgb(1, 0, 0, 1/4))
ecdfPlot(
  asset[asset > quantile(asset, probs = .95, na.rm = T)] - 
         quantile(asset, probs = .95, na.rm = T), add = T, ecdf.col = rgb(0, 1, 0, 1/4))

##———ОБЛИГАЦИИ————

ass_num = 4

asset = all_data[, 4] %>% drop_na() %>% unlist() %>% as.numeric()

asset.log = log(asset - min(asset))
dens = density(asset.log, kernel = "optcosine", adjust = 2)
cdff = CDF(dens)
invcdf <- function(q) {
  uniroot(function(x) {cdff(x) - q}, range(asset.log), lower = -100)$root
}
icdf.bonds = to_vec(for (i in x[,ass_num]) exp(as.numeric(try(invcdf(i), silent = TRUE))) + min(asset))


#dens = density(asset, kernel = "optcosine", adjust = 3.5)
#cdff = CDF(dens)
#invcdf <- function(q) {
#  uniroot(function(x) {cdff(x) - q}, range(asset), lower = -100)$root
#}
#icdf.bonds = to_vec(for (i in x[,ass_num]) as.numeric(try(invcdf(i), silent = TRUE)))
p = 0.05
n = as.integer(asset %>% length() * p)

pareto.lower = gpdFit(-asset, nextremes = n, method = "mle")
pareto.upper = gpdFit(asset, nextremes = n, method = "mle")

upper.loc = exp(as.numeric(try(invcdf(1-p), silent = TRUE))) + min(asset)
lower.loc = -(exp(as.numeric(try(invcdf(p), silent = TRUE))) + min(asset))

#upper.loc = as.numeric(try(invcdf(1-p), silent = TRUE))
#lower.loc = -as.numeric(try(invcdf(p), silent = TRUE))

bonds = to_vec(
  for (i in x[,ass_num])
    if (i >= (1 - p)) as.numeric(try(qgpd((i - (1 - p))/(max(x[,ass_num][x[,ass_num] > (1 - p)]) - 
                                                           min(x[,ass_num][x[,ass_num] > (1 - p)])), 
                                          loc = upper.loc, 
                                          scale = pareto.upper$par.ests[1], 
                                          shape = pareto.upper$par.ests[2]), silent = TRUE))
    else if (i <= p) as.numeric(try(-qgpd(i/(max(x[,ass_num][x[,ass_num] <= p]) - 
                                               min(x[,ass_num][x[,ass_num] <= p])), 
                                            loc = lower.loc, 
                                            scale = pareto.lower$par.ests[1], 
                                            shape = pareto.lower$par.ests[2]), silent = TRUE))
    else exp(as.numeric(try(invcdf(i), silent = TRUE))) + min(asset)
    #else as.numeric(try(invcdf(i), silent = TRUE))
  )

hist(asset, breaks = 100, col = rgb(0, 0, 1, 1/4), freq = F)
hist(icdf.bonds, breaks = 100, add = T, col = rgb(1, 0, 0, 1/4), freq = F)
hist(bonds, 
     breaks = 300,
     add = T, 
     col = rgb(0, 1, 0, 1/4), freq = F)

y1 <- year.ret(bonds, 1)
y1 %>% median(na.rm = T)
## [1] 0.09820706
y1 %>% quantile(na.rm = T, probs = c(0.001, 0.01, .1, 0.9, 0.99, .999))
##        0.1%          1%         10%         90%         99%       99.9% 
## -0.01981055  0.01231070  0.05149279  0.14607927  0.19034115  0.22252871
y1.log <- year.ret(icdf.bonds, 1)
y1.log %>% median(na.rm = T)
## [1] 0.0981027
y1.log %>% quantile(na.rm = T, probs = c(0.001, 0.01, .1, 0.9, 0.99, .999))
##         0.1%           1%          10%          90%          99%        99.9% 
## -0.009849212  0.016656062  0.054501748  0.145013423  0.188142549  0.217129775
y1.h <- year.ret(asset, 1)
y1.h %>% median(na.rm = T)
## [1] 0.09591952
y1.h %>% quantile(na.rm = T, probs = c(0.001, 0.01, .1, 0.9, 0.99, .999))
##        0.1%          1%         10%         90%         99%       99.9% 
## -0.08930894 -0.06569672  0.03337894  0.17933939  0.31541295  0.34104698
hist(y1.h, 
     freq = F, breaks = 100, col = rgb(1, 0, 0, 1/4))
hist(y1, freq = F, breaks = 150, add = T, col = rgb(0, 1, 0, 1/4))

hist(y1.h, 
     freq = F, breaks = 100, col = rgb(1, 0, 0, 1/4))
hist(y1.log, freq = F, breaks = 100, add = T, col = rgb(0, 1, 0, 1/4))

ecdfPlot(
  bonds[bonds > quantile(bonds, probs = .95, na.rm = T)] - 
         quantile(bonds, probs = .95, na.rm = T), ecdf.col = rgb(1, 0, 0, 1/4))
ecdfPlot(
  asset[asset > quantile(asset, probs = .95, na.rm = T)] - 
         quantile(asset, probs = .95, na.rm = T), add = T, ecdf.col = rgb(0, 1, 0, 1/4))

ecdfPlot(
  icdf.bonds[icdf.bonds > quantile(icdf.bonds, probs = .95, na.rm = T)] - 
         quantile(icdf.bonds, probs = .95, na.rm = T), ecdf.col = rgb(1, 0, 0, 1/4))
ecdfPlot(
  asset[asset > quantile(asset, probs = .95, na.rm = T)] - 
         quantile(asset, probs = .95, na.rm = T), add = T, ecdf.col = rgb(0, 1, 0, 1/4))

##———ДОЛЛАР————

ass_num  = 2
asset = unlist(all_data[, ass_num] %>% drop_na()) %>% as.numeric
asset.log = log(asset - min(asset))
dens = density(asset.log, kernel = "optcosine", adjust = 2)
cdff = CDF(dens)
invcdf <- function(q) {
  uniroot(function(x) {cdff(x) - q}, range(asset.log), lower = -100)$root
}

icdf.usd = to_vec(for (i in x[,2]) exp(as.numeric(try(invcdf(i), silent = TRUE))) + min(asset))

hist(icdf.usd, breaks = 100, freq = F, col = rgb(1, 0, 0, 1/4))
hist(asset, breaks = 100, freq = F, add = T, col = rgb(0, 1, 0, 1/4))

p = 0.05

pareto.lower = gpdFit(-asset, threshold = p, method = "mle")
pareto.upper = gpdFit(asset, threshold = p, method = "mle")

upper.loc = exp(as.numeric(try(invcdf(1-p), silent = TRUE))) + min(asset)
lower.loc = -(exp(as.numeric(try(invcdf(p), silent = TRUE))) + min(asset))

usd = to_vec(
  for (i in x[,ass_num]) 
    if (i >= (1 - p)) as.numeric(try(qgpd((i - (1 - p))/(max(x[,ass_num][x[,ass_num] > (1 - p)]) - 
                                                           min(x[,ass_num][x[,ass_num] > (1 - p)])), 
                                          loc = upper.loc, 
                                          scale = pareto.upper$par.ests[1], 
                                          shape = pareto.upper$par.ests[2]), silent = TRUE))
    else if (i <= p) as.numeric(try(-qgpd(i/(max(x[,ass_num][x[,ass_num] <= p]) - 
                                               min(x[,ass_num][x[,ass_num] <= p])), 
                                            loc = lower.loc, 
                                            scale = pareto.lower$par.ests[1], 
                                            shape = pareto.lower$par.ests[2]), silent = TRUE))
    else exp(as.numeric(try(invcdf(i), silent = TRUE))) + min(asset)
  )

hist(asset, breaks = 200, col = rgb(0, 0, 1, 1/4), freq = F)
hist(icdf.usd, breaks = 100, add = T, col = rgb(1, 0, 0, 1/4), freq = F)
hist(usd, 
     breaks = 300,
     add = T, 
     col = rgb(0, 1, 0, 1/4), freq = F)

y1 <- year.ret(usd, 1)
y1 %>% median(na.rm = T)
## [1] 0.09413204
y1 %>% quantile(na.rm = T, probs = c(0.001, 0.01, .1, 0.9, 0.99, .999))
##       0.1%         1%        10%        90%        99%      99.9% 
## -0.4003052 -0.2898254 -0.1386057  0.3772739  0.6469161  0.9029923
y1.log <- year.ret(icdf.usd, 1)
y1.log %>% median(na.rm = T)
## [1] 0.09948863
y1.log %>% quantile(na.rm = T, probs = c(0.001, 0.01, .1, 0.9, 0.99, .999))
##       0.1%         1%        10%        90%        99%      99.9% 
## -0.3829171 -0.2789847 -0.1272547  0.3786651  0.6398721  0.8932549
y1.h <- year.ret(asset, 1)
y1.h %>% median(na.rm = T)
## [1] 0.1029372
y1.h %>% quantile(na.rm = T, probs = c(0.001, 0.01, .1, 0.9, 0.99, .999))
##        0.1%          1%         10%         90%         99%       99.9% 
## -0.27856019 -0.24939644 -0.07849765  0.32471088  0.49785623  0.59729669
hist(y1.h, 
     freq = F, breaks = 100, col = rgb(1, 0, 0, 1/4))
hist(y1, freq = F, breaks = 100, add = T, col = rgb(0, 1, 0, 1/4))

hist(y1.h, 
     freq = F, breaks = 100, col = rgb(1, 0, 0, 1/4))
hist(y1.log, freq = F, breaks = 100, add = T, col = rgb(0, 1, 0, 1/4))

ecdfPlot(
  usd[usd > quantile(usd, probs = .95, na.rm = T)] - 
         quantile(usd, probs = .95, na.rm = T), ecdf.col = rgb(1, 0, 0, 1/4))
ecdfPlot(
  asset[asset > quantile(asset, probs = .95, na.rm = T)] - 
         quantile(asset, probs = .95, na.rm = T), add = T, ecdf.col = rgb(0, 1, 0, 1/4))

ecdfPlot(
  icdf.usd[icdf.usd > quantile(icdf.usd, probs = .95, na.rm = T)] - 
         quantile(icdf.usd, probs = .95, na.rm = T), ecdf.col = rgb(1, 0, 0, 1/4))
ecdfPlot(
  asset[asset > quantile(asset, probs = .95, na.rm = T)] - 
         quantile(asset, probs = .95, na.rm = T), add = T, ecdf.col = rgb(0, 1, 0, 1/4))

##———ЗОЛОТО ДОЛЛ————

ass_num  = 1
asset = unlist(all_data[, ass_num] %>% drop_na()) %>% as.numeric
asset.log = log(asset - min(asset))
dens = density(asset, kernel = "optcosine", adjust = 2)
cdff = CDF(dens)
invcdf <- function(q) {
  uniroot(function(x) {cdff(x) - q}, range(asset), lower = -100)$root
}

#icdf.gold = to_vec(for (i in x[,ass_num]) exp(as.numeric(try(invcdf(i), silent = TRUE))) + min(asset))
icdf.gold = to_vec(for (i in x[,ass_num]) as.numeric(try(invcdf(i), silent = TRUE)))
p = 0.05

pareto.lower = gpdFit(-asset, threshold = p, method = "mle")
pareto.upper = gpdFit(asset, threshold = p, method = "mle")

#upper.loc = exp(as.numeric(try(invcdf(1-p), silent = TRUE))) + min(asset)
#lower.loc = -(exp(as.numeric(try(invcdf(p), silent = TRUE))) + min(asset))

upper.loc = as.numeric(try(invcdf(1-p), silent = TRUE))
lower.loc = -as.numeric(try(invcdf(p), silent = TRUE))

gold = to_vec(
  for (i in x[,ass_num]) 
    if (i >= (1 - p)) as.numeric(try(qgpd((i - (1 - p))/(max(x[,ass_num][x[,ass_num] > (1 - p)]) - 
                                                           min(x[,ass_num][x[,ass_num] > (1 - p)])), 
                                          loc = upper.loc, 
                                          scale = pareto.upper$par.ests[1], 
                                          shape = pareto.upper$par.ests[2]), silent = TRUE))
    else if (i <= p) as.numeric(try(-qgpd(i/(max(x[,ass_num][x[,ass_num] <= p]) - 
                                               min(x[,ass_num][x[,ass_num] <= p])), 
                                            loc = lower.loc, 
                                            scale = pareto.lower$par.ests[1], 
                                            shape = pareto.lower$par.ests[2]), silent = TRUE))
    #else exp(as.numeric(try(invcdf(i), silent = TRUE))) + min(asset)
    else as.numeric(try(invcdf(i), silent = TRUE))
  )

hist(asset, breaks = 200, col = rgb(0, 0, 1, 1/4), freq = F)
hist(icdf.gold, breaks = 100, add = T, col = rgb(1, 0, 0, 1/4), freq = F)
hist(gold, 
     breaks = 300,
     add = T, 
     col = rgb(0, 1, 0, 1/4), freq = F)

y1 <- year.ret(gold, 1)
y1 %>% median(na.rm = T)
## [1] 0.1317747
y1 %>% quantile(na.rm = T, probs = c(0.001, 0.01, .1, 0.9, 0.99, .999))
##       0.1%         1%        10%        90%        99%      99.9% 
## -0.4826304 -0.3517421 -0.1637771  0.5206318  0.9491465  1.3020972
y1.log <- year.ret(icdf.gold, 1)
y1.log %>% median(na.rm = T)
## [1] 0.1177174
y1.log %>% quantile(na.rm = T, probs = c(0.001, 0.01, .1, 0.9, 0.99, .999))
##       0.1%         1%        10%        90%        99%      99.9% 
## -0.4420908 -0.3507456 -0.1748621  0.5037450  0.9006297  1.2662826
y1.h <- year.ret(asset, 1)
y1.h %>% median(na.rm = T)
## [1] 0.1634978
y1.h %>% quantile(na.rm = T, probs = c(0.001, 0.01, .1, 0.9, 0.99, .999))
##       0.1%         1%        10%        90%        99%      99.9% 
## -0.2830548 -0.2369268 -0.1150337  0.4846184  0.7980326  0.9731775
hist(y1, 
     freq = F, breaks = 100, col = rgb(1, 0, 0, 1/4))
hist(y1.h, freq = F, breaks = 100, add = T, col = rgb(0, 1, 0, 1/4))

hist(y1.log, 
     freq = F, breaks = 100, col = rgb(1, 0, 0, 1/4))
hist(y1.h, freq = F, breaks = 100, add = T, col = rgb(0, 1, 0, 1/4))

ecdfPlot(
  gold[gold > quantile(gold, probs = .95, na.rm = T)] - 
         quantile(gold, probs = .95, na.rm = T), ecdf.col = rgb(1, 0, 0, 1/4))
ecdfPlot(
  asset[asset > quantile(asset, probs = .95, na.rm = T)] - 
         quantile(asset, probs = .95, na.rm = T), add = T, ecdf.col = rgb(0, 1, 0, 1/4))

ecdfPlot(
  icdf.gold[icdf.gold > quantile(icdf.gold, probs = .95, na.rm = T)] - 
         quantile(icdf.gold, probs = .95, na.rm = T), ecdf.col = rgb(1, 0, 0, 1/4))
ecdfPlot(
  asset[asset > quantile(asset, probs = .95, na.rm = T)] - 
         quantile(asset, probs = .95, na.rm = T), add = T, ecdf.col = rgb(0, 1, 0, 1/4))

##———ЗОЛОТО————

ass_num  = 5
asset = unlist(all_data[, ass_num] %>% drop_na()) %>% as.numeric
asset.log = log(asset - min(asset))
gold.rus = (1 + gold) * (1 + usd) - 1
icdf.gold.rus = (1 + icdf.gold) * (1 + icdf.usd) - 1
gold.rus2 = (1 + icdf.gold) * (1 + usd) - 1

hist(gold.rus, freq = F, breaks = 100, col = rgb(1, 0, 0, 1/4))
hist(icdf.gold.rus, freq = F, add = T, breaks = 100, col = rgb(0, 0, 1, 1/4))
hist(asset, freq = F, add = T, breaks = 100, col = rgb(0, 1, 0, 1/4))

y1 <- year.ret(gold.rus, 1)
y1 %>% median(na.rm = T)
## [1] 0.2343543
y1 %>% quantile(na.rm = T, probs = c(0.001, 0.01, .1, 0.9, 0.99, .999))
##       0.1%         1%        10%        90%        99%      99.9% 
## -0.4971333 -0.3631454 -0.1333392  0.7425511  1.3716187  2.1124551
y1.2 <- year.ret(gold.rus2, 1)
y1.2 %>% median(na.rm = T)
## [1] 0.2198772
y1.2 %>% quantile(na.rm = T, probs = c(0.001, 0.01, .1, 0.9, 0.99, .999))
##       0.1%         1%        10%        90%        99%      99.9% 
## -0.4712816 -0.3549255 -0.1415999  0.7120560  1.3243026  2.0714529
y1.log <- year.ret(icdf.gold.rus, 1)
y1.log %>% median(na.rm = T)
## [1] 0.2289444
y1.log %>% quantile(na.rm = T, probs = c(0.001, 0.01, .1, 0.9, 0.99, .999))
##       0.1%         1%        10%        90%        99%      99.9% 
## -0.4665207 -0.3430369 -0.1302717  0.7211038  1.3245494  2.0630525
y1.h <- year.ret(asset, 1)
y1.h %>% median(na.rm = T)
## [1] 0.019807
y1.h %>% quantile(na.rm = T, probs = c(0.001, 0.01, .1, 0.9, 0.99, .999))
##        0.1%          1%         10%         90%         99%       99.9% 
## -0.26386767 -0.21837821 -0.09742842  0.35514442  0.80059683  0.93646183
hist(y1, 
     freq = F, breaks = 100, col = rgb(1, 0, 0, 1/4))
hist(y1.h, freq = F, breaks = 30, add = T, col = rgb(0, 1, 0, 1/4))

hist(y1.2, 
     freq = F, breaks = 100, col = rgb(1, 0, 0, 1/4))
hist(y1.h, freq = F, breaks = 30, add = T, col = rgb(0, 1, 0, 1/4))

hist(y1.log, 
     freq = F, breaks = 100, col = rgb(1, 0, 0, 1/4))
hist(y1.h, freq = F, breaks = 30, add = T, col = rgb(0, 1, 0, 1/4))

ecdfPlot(
  gold.rus2[gold.rus2 > quantile(gold.rus2, probs = .95, na.rm = T)] - 
         quantile(gold.rus2, probs = .95, na.rm = T), ecdf.col = rgb(0, 0, 1, 1/2))
ecdfPlot(
  asset[asset > quantile(asset, probs = .95, na.rm = T)] - 
         quantile(asset, probs = .95, na.rm = T), add = T, ecdf.col = rgb(0, 1, 0, 1/2))

ecdfPlot(
  gold.rus[gold.rus > quantile(gold.rus, probs = .95, na.rm = T)] - 
         quantile(gold.rus, probs = .95, na.rm = T), ecdf.col = rgb(0, 0, 1, 1/2))
ecdfPlot(
  asset[asset > quantile(asset, probs = .95, na.rm = T)] - 
         quantile(asset, probs = .95, na.rm = T), add = T, ecdf.col = rgb(0, 1, 0, 1/2))

ecdfPlot(
  icdf.gold.rus[icdf.gold.rus > quantile(icdf.gold.rus, probs = .95, na.rm = T)] - 
         quantile(icdf.gold.rus, probs = .95, na.rm = T), ecdf.col = rgb(1, 0, 0, 1/4))
ecdfPlot(
  asset[asset > quantile(asset, probs = .95, na.rm = T)] - 
         quantile(asset, probs = .95, na.rm = T), add = T, ecdf.col = rgb(0, 1, 0, 1/2))

##———ИНФЛЯЦИЯ————

ass_num  = 5
asset = (unlist(all_data[1074:1277, ass_num] %>% drop_na()) %>% as.numeric)

asset <- asset[!is.na(asset)]
asset.log = log(asset - min(asset))
if (ass_num == 6){
  ass_num = 5
}

dens = density(asset, kernel = "optcosine", adjust = 2)
cdff = CDF(dens)
invcdf <- function(q) {
  uniroot(function(x) {cdff(x) - q}, range(asset), lower = -100)$root
}

#icdf.infl = to_vec(for (i in x[,ass_num]) exp(as.numeric(try(invcdf(i), silent = TRUE))) + min(asset))
icdf.infl = to_vec(for (i in x[,ass_num]) as.numeric(try(invcdf(i), silent = TRUE)))
p = 0.05
n = as.integer(asset %>% length() * p)

pareto.lower = gpdFit(-asset, nextremes = 30, method = "mle")
pareto.upper = gpdFit(asset, nextremes = 20, method = "mle")

#upper.loc = exp(as.numeric(try(invcdf(1-p), silent = TRUE))) + min(asset)
#lower.loc = -(exp(as.numeric(try(invcdf(p), silent = TRUE))) + min(asset))

upper.loc = as.numeric(try(invcdf(1-p), silent = TRUE))
lower.loc = -as.numeric(try(invcdf(p), silent = TRUE))

infl = to_vec(
  for (i in x[,ass_num])
    if (i >= (1 - p)) as.numeric(try(qgpd((i - (1 - p))/(max(x[,ass_num][x[,ass_num] > (1 - p)]) - 
                                                           min(x[,ass_num][x[,ass_num] > (1 - p)])), 
                                          loc = upper.loc, 
                                          scale = pareto.upper$par.ests[1], 
                                          shape = pareto.upper$par.ests[2]), silent = TRUE))
    else if (i <= p) as.numeric(try(-qgpd(i/(max(x[,ass_num][x[,ass_num] <= p]) - 
                                               min(x[,ass_num][x[,ass_num] <= p])), 
                                            loc = lower.loc, 
                                            scale = pareto.lower$par.ests[1], 
                                            shape = pareto.lower$par.ests[2]), silent = TRUE))
    #else exp(as.numeric(try(invcdf(i), silent = TRUE))) + min(asset)
    else as.numeric(try(invcdf(i), silent = TRUE))
  )

hist(asset, breaks = 200, col = rgb(0, 0, 1, 1/4), freq = F)
hist(icdf.infl, breaks = 100, 
     add = T, 
     col = rgb(1, 0, 0, 1/4), freq = F)
hist(infl, 
     breaks = 300,
     add = T, 
     col = rgb(0, 1, 0, 1/4), freq = F)

-as.numeric(try(invcdf(p), silent = TRUE))
## [1] 0.05386012
hist(asset[asset <= -lower.loc], freq = F, breaks = 100, col = rgb(0, 0, 1, 1/4))

hist(asset[asset > upper.loc], freq = F, breaks = 100, col = rgb(0, 0, 1, 1/4))

hist(rgpd(1000, 
          loc = upper.loc, 
          scale = pareto.upper$par.ests[1], 
          shape = pareto.upper$par.ests[2]), breaks = 100)

hist(-rgpd(1000, 
          loc = lower.loc, 
          scale = pareto.lower$par.ests[1]*10000, 
          shape = pareto.lower$par.ests[2]/2), breaks = 100)

hist(infl, freq = F, col = rgb(1, 0, 0, 1/4), breaks = 100)
hist(asset, freq = F, add = T, col = rgb(0, 1, 0, 1/4), breaks = 100)

hist(icdf.infl, freq = F, col = rgb(0, 0, 1, 1/4), breaks = 100, add = T)
lines(dens)

y1 <- year.ret(infl[abs(infl) <  0.01], 1)
y1 %>% mean(na.rm = T)
## [1] -0.004045221
y1 %>% median(na.rm = T)
## [1] -0.005037082
y1 %>% quantile(na.rm = T, probs = c(0.001, 0.01, .1, 0.9, 0.99, .999))
##        0.1%          1%         10%         90%         99%       99.9% 
## -0.11383358 -0.09360653 -0.05285529  0.04727358  0.08924713  0.12746680
y1.log <- year.ret(icdf.infl, 1)
y1.log %>% mean(na.rm = T)
## [1] 0.009530531
y1.log %>% median(na.rm = T)
## [1] -0.02202364
y1.log %>% quantile(na.rm = T, probs = c(0.001, 0.01, .1, 0.9, 0.99, .999))
##       0.1%         1%        10%        90%        99%      99.9% 
## -0.5444674 -0.4480876 -0.2828437  0.3454748  0.7357613  1.0024144
y1.h <- year.ret(asset, 1)
y1.h %>% mean(na.rm = T)
## [1] 0.108378
y1.h %>% median(na.rm = T)
## [1] 0.03463847
y1.h %>% quantile(na.rm = T, probs = c(0.001, 0.01, .1, 0.9, 0.99, .999))
##       0.1%         1%        10%        90%        99%      99.9% 
## -0.2624147 -0.2405996 -0.1591392  0.5112031  0.6258575  0.7260445
hist(y1, 
     freq = F, breaks = 100, col = rgb(1, 0, 0, 1/4))
hist(y1.h, freq = F, breaks = 100, add = T, col = rgb(0, 1, 0, 1/4))

hist(y1.log,
     freq = F, breaks = 100, col = rgb(1, 0, 0, 1/4))
hist(y1.h, freq = F, breaks = 100,
     add = T,
     col = rgb(0, 1, 0, 1/4))

ecdfPlot(
  infl[infl > quantile(infl, probs = .95, na.rm = T)] -
         quantile(infl, probs = .95, na.rm = T), ecdf.col = rgb(0, 0, 1, 1/2))
## Warning in is.not.finite.warning(x): There were 2 nonfinite values in x : 2
## NA's
## Warning in ecdfPlot(infl[infl > quantile(infl, probs = 0.95, na.rm = T)] - : 2
## observations with NA/NaN/Inf in 'x' removed.
ecdfPlot(
  asset[asset > quantile(asset, probs = .95, na.rm = T)] -
         quantile(asset, probs = .95, na.rm = T), add = T, 
  ecdf.col = rgb(0, 1, 0, 1/2))

ecdfPlot(
  icdf.infl[icdf.infl > quantile(icdf.infl, probs = .95, na.rm = T)] - 
         quantile(icdf.infl, probs = .95, na.rm = T), ecdf.col = rgb(0, 0, 1, 1/2))
## Warning in is.not.finite.warning(x): There were 262 nonfinite values in x : 262
## NA's
## Warning in ecdfPlot(icdf.infl[icdf.infl > quantile(icdf.infl, probs = 0.95, :
## 262 observations with NA/NaN/Inf in 'x' removed.
ecdfPlot(
  asset[asset > quantile(asset, probs = .95, na.rm = T)] - 
         quantile(asset, probs = .95, na.rm = T), add = T, ecdf.col = rgb(0, 1, 0, 1/2))

ecdfPlot(
  infl[infl > quantile(infl, probs = .95, na.rm = T)] - 
         quantile(infl, probs = .95, na.rm = T), ecdf.col = rgb(0, 0, 1, 1/2))
ecdfPlot(
  asset[asset > quantile(asset, probs = .95, na.rm = T)] - 
         quantile(asset, probs = .95, na.rm = T), add = T, ecdf.col = rgb(0, 1, 0, 1/2))

ecdfPlot(
  icdf.infl[icdf.infl > quantile(icdf.infl, probs = .95, na.rm = T)] - 
         quantile(icdf.infl, probs = .95, na.rm = T), ecdf.col = rgb(0, 0, 1, 1/2))
ecdfPlot(
  asset[asset > quantile(asset, probs = .95, na.rm = T)] - 
         quantile(asset, probs = .95, na.rm = T), add = T, ecdf.col = rgb(0, 1, 0, 1/2))

result <- data.frame(gold, usd, equity + 0.0007, bonds, gold.rus2, infl) %>% drop_na()
colnames(result) <- c("gold", "usd", "equity", "bonds", "gold.rus", "infl")
result.real <- (result[,3:5] + 1)/(result[,6] + 1) - 1

all_data %>% drop_na() %>% cor()
##                  Gold_return Gold_USD_return RusEquity_return RusBonds_return
## Gold_return       1.00000000      0.67233415       0.04379932     -0.05742696
## Gold_USD_return   0.67233415      1.00000000       0.19762326      0.09615949
## RusEquity_return  0.04379932      0.19762326       1.00000000      0.26418438
## RusBonds_return  -0.05742696      0.09615949       0.26418438      1.00000000
## USD_return        0.64715627     -0.12886925      -0.14504625     -0.18072901
##                  USD_return
## Gold_return       0.6471563
## Gold_USD_return  -0.1288692
## RusEquity_return -0.1450462
## RusBonds_return  -0.1807290
## USD_return        1.0000000
result.real %>% cor()
##             equity     bonds  gold.rus
## equity   1.0000000 0.7679124 0.5326608
## bonds    0.7679124 1.0000000 0.6650153
## gold.rus 0.5326608 0.6650153 1.0000000
x %>% cor()
##             [,1]        [,2]        [,3]        [,4]        [,5]
## [1,]  1.00000000 -0.16556989  0.12064037  0.05042423  0.09866808
## [2,] -0.16556989  1.00000000 -0.09038438 -0.19648459 -0.14238535
## [3,]  0.12064037 -0.09038438  1.00000000  0.25268633 -0.02873913
## [4,]  0.05042423 -0.19648459  0.25268633  1.00000000 -0.07326493
## [5,]  0.09866808 -0.14238535 -0.02873913 -0.07326493  1.00000000
hist(result$infl, breaks = 200)

y1.eq <- year.ret(result.real[,1], 1)
y1.eq %>% median(na.rm = T)
## [1] 0.1776386
y1.eq %>% mean(na.rm = T)
## [1] 0.245799
y1.bnd <- year.ret(result.real[,2], 1)
y1.bnd %>% median(na.rm = T)
## [1] 0.09760559
y1.bnd %>% mean(na.rm = T)
## [1] 0.1328201
y1.gld <- year.ret(result.real[,3], 1)
y1.gld %>% median(na.rm = T)
## [1] 0.2288711
y1.gld %>% mean(na.rm = T)
## [1] 0.3020637
write.csv(data.frame(gold, usd, equity + 0.0006, bonds, gold.rus) %>% drop_na(),
          "modeled_R.csv", 
          row.names = FALSE)

write.csv(result.real %>% drop_na(),
          "modeled_R_real2.csv", 
          row.names = FALSE)

##ContourPlots

df <- data.frame(gold = x[, 1],
                 usd = x[, 2],
                 equity = x[, 3],
                 bonds = x[, 4], 
                 infl = x[, 5])

df.mod <- data.frame(gold,
                     usd,
                     equity,
                     bonds,
                     infl) %>% 
  drop_na()

df.mod.icdf <- data.frame(gold = icdf.gold,
                     usd = icdf.usd,
                     equity = icdf.equity,
                     bonds = icdf.bonds,
                     gold.rus = icdf.infl) %>% 
  drop_na()
df.mod.real <- result %>% cbind(result.real)
colnames(df.mod.real)[7:9] <- c("equity.real", "bonds.real", "gold.rus.real") 

df.mod.real <- as.data.frame(apply(
  df.mod.real[,c("equity", "bonds", "gold.rus", "equity.real", "bonds.real", "gold.rus.real")], 
  2, 
  year.ret))
ggplot(df.mod.real) +
  geom_histogram(aes(x = equity), fill="#69b3a2", alpha = .5, bins = 100) +
  geom_histogram(aes(x = equity.real), fill = '#e95c3d', alpha = .5, bins = 100) +
  geom_vline(aes(xintercept = median(equity))) +
  geom_vline(aes(xintercept = median(equity.real)))

ggplot(df.mod.real) +
  geom_histogram(aes(x = bonds), fill="#69b3a2", alpha = .5, bins = 100) +
  geom_histogram(aes(x = bonds.real), fill = '#e95c3d', alpha = .5, bins = 100) +
  geom_vline(aes(xintercept = median(bonds))) +
  geom_vline(aes(xintercept = median(bonds.real)))

ggplot(df.mod.real) +
  geom_histogram(aes(x = gold.rus), fill="#69b3a2", alpha = .5, bins = 100) +
  geom_histogram(aes(x = gold.rus.real), fill = '#e95c3d', alpha = .5, bins = 100) +
  geom_vline(aes(xintercept = median(gold.rus))) +
  geom_vline(aes(xintercept = median(gold.rus.real)))

df.mod.real$equity %>% median()
## [1] 0.1817051
df.mod.real$bonds %>% median()
## [1] 0.09797723
df.mod.real$gold.rus %>% median()
## [1] 0.219557
df.mod.real$equity.real %>% median()
## [1] 0.1776386
df.mod.real$bonds.real %>% median()
## [1] 0.09760559
df.mod.real$gold.rus.real %>% median()
## [1] 0.2288711
ggplot(df.mod.real, aes(x = equity, y = equity.real)) +
  geom_density_2d(aes(color = after_stat(level))) +
  scale_color_viridis_c()

ggplot(df.mod.real, aes(x = bonds, y = bonds.real)) +
  geom_density_2d(aes(color = after_stat(level))) +
  scale_color_viridis_c()

ggplot(df.mod.real, aes(x = gold.rus, y = gold.rus.real)) +
  geom_density_2d(aes(color = after_stat(level))) +
  scale_color_viridis_c()

ggplot(df.mod.icdf, aes(x = gold, y = usd)) +
  geom_density_2d(aes(color = ..level..)) +
  scale_color_viridis_c()
## Warning: The dot-dot notation (`..level..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(level)` instead.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

ggplot(df.mod.icdf, aes(x = equity, y = bonds)) +
  geom_density_2d(aes(color = ..level..)) +
  scale_color_viridis_c()

ggplot(df.mod.icdf, aes(x = equity, y = gold)) +
  geom_density_2d(aes(color = ..level..)) +
  scale_color_viridis_c()

ggplot(df.mod.icdf, aes(x = equity, y = usd)) +
  geom_density_2d(aes(color = ..level..)) +
  scale_color_viridis_c()

ggplot(df, aes(x = gold, y = usd)) +
  geom_density_2d_filled(alpha = 0.75) +
  scale_color_viridis_c()

ggplot(df, aes(x = equity, y = bonds)) +
  geom_density_2d_filled(alpha = 0.75) +
  scale_color_viridis_c()

ggplot(df, aes(x = equity, y = gold)) +
  geom_density_2d_filled(alpha = 0.75) +
  scale_color_viridis_c()

ggplot(df, aes(x = equity, y = usd)) +
  geom_density_2d_filled(alpha = 0.75) +
  scale_color_viridis_c()