getwd()
## [1] "F:/THESIS/Tanya/R-Folder/data"
library(readxl)
## Warning: package 'readxl' was built under R version 3.4.4
goy <- read_excel("goy.xls")
## readxl works best with a newer version of the tibble package.
## You currently have tibble v1.4.2.
## Falling back to column name repair from tibble <= v1.4.2.
## Message displays once per session.
goy$observation_date<- as.Date(goy$observation_date)
str(goy)
## Classes 'tbl_df', 'tbl' and 'data.frame':    879 obs. of  4 variables:
##  $ observation_date: Date, format: "1946-01-01" "1946-02-01" ...
##  $ gold            : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ oil             : num  1.17 1.17 1.17 1.27 1.27 1.27 1.27 1.52 1.52 1.52 ...
##  $ yen             : num  NA NA NA NA NA NA NA NA NA NA ...
library(FRAPO)
## Warning: package 'FRAPO' was built under R version 3.4.4
## Loading required package: cccp
## Warning: package 'cccp' was built under R version 3.4.4
## Loading required package: Rglpk
## Warning: package 'Rglpk' was built under R version 3.4.4
## Loading required package: slam
## Warning: package 'slam' was built under R version 3.4.4
## Using the GLPK callable library version 4.47
## Loading required package: timeSeries
## Warning: package 'timeSeries' was built under R version 3.4.4
## Loading required package: timeDate
## Warning: package 'timeDate' was built under R version 3.4.4
## Financial Risk Modelling and Portfolio Optimisation with R (version 0.4-1)
library(timeSeries)
library(QRM)
## Warning: package 'QRM' was built under R version 3.4.4
## Loading required package: gsl
## Warning: package 'gsl' was built under R version 3.4.4
## Loading required package: Matrix
## Loading required package: mvtnorm
## Warning: package 'mvtnorm' was built under R version 3.4.4
## Loading required package: numDeriv
## Warning: package 'numDeriv' was built under R version 3.4.4
## 
## Attaching package: 'QRM'
## The following object is masked from 'package:base':
## 
##     lbeta
library(fGarch)
## Warning: package 'fGarch' was built under R version 3.4.4
## Loading required package: fBasics
## Warning: package 'fBasics' was built under R version 3.4.4
goycc<- na.omit(goy)

date<- goycc$observation_date
str(date)
##  Date[1:578], format: "1971-01-01" "1971-02-01" "1971-03-01" "1971-04-01" "1971-05-01" ...
date<- as.character(date)
str(date)
##  chr [1:578] "1971-01-01" "1971-02-01" "1971-03-01" "1971-04-01" ...
goydata<- subset(goycc[2:4])
goyts<- timeSeries(goydata,charvec = date, title = "goyts")
goydatats<-timeSeries(goydata)
goyloss<-as.data.frame(na.omit(-1.0*diff(log(goydatats))*100.0))
head(goyloss)
##         gold oil          yen
## 1 -2.2159727   0 0.1327622424
## 2 -0.4042497   0 0.0073572041
## 3 -0.3291838   0 0.0043394012
## 4 -3.7529485   0 0.0252286411
## 5  0.9621424   0 0.0003306595
## 6 -2.0351284   0 0.0021075260
acf(goyts$gold)

acf(goyts$oil)

acf(goyts$yen)

pacf(goyts$gold)

pacf(goyts$oil)

pacf(goyts$yen)

arimagold<- arima(goyts$gold)
arimaoil<- arima(goyts$oil)
arimayen<- arima(goyts$yen)
ggold<- garchFit(formula = ~garch(1,1), data = goydatats$gold, cond.dist = "std", trace = FALSE)
## Warning in sqrt(diag(fit$cvar)): NaNs produced
coef(ggold)
##          mu       omega      alpha1       beta1       shape 
## 382.6841257  23.4439867   0.8340353   0.3549226  10.0000000
plot(sqrt(252) * ggold@sigma.t, type="l")

goil<- garchFit(formula = ~garch(1,1), data = goydatats$oil, cond.dist = "std", trace = FALSE)
coef(goil)
##           mu        omega       alpha1        beta1        shape 
## 1.485167e+01 7.695921e-04 1.000000e+00 2.354893e-01 9.689374e+00
plot(sqrt(252) * goil@sigma.t, type="l")

gyen<- garchFit(formula = ~garch(1,1), data = goydatats$yen, cond.dist = "std", trace = FALSE)
## Warning in sqrt(diag(fit$cvar)): NaNs produced
coef(gyen)
##          mu       omega      alpha1       beta1       shape 
## 110.8960267   4.1177145   1.0000000   0.1242049  10.0000000
plot(sqrt(252) * gyen@sigma.t, type="l")

goldstd<- stdev(goydata$gold)
goldstd
## [1] 428.8314
oilstd<- stdev(goydata$oil)
oilstd
## [1] 27.74152
yenstd<- stdev(goydata$yen)
yenstd
## [1] 73.12937
ESgarchg <- function(y, p = 0.99){
  gfit <- garchFit(formula = ~garch(1, 1), data = y$gold,
                   cond.dist = "std", trace = FALSE)
  sigma <-  predict(gfit, n.ahead = 1)[3]
  df <- coef(gfit)["shape"]
  ES <- sigma * (dt(qt(p, df), df)/(1 - p)) *
    ((df + (qt(p, df))^2)/(df - 1))
  return(ES)
}
ESgarchg(goydatats,p=0.95)
## Warning in sqrt(diag(fit$cvar)): NaNs produced
##   standardDeviation
## 1          2521.724
ESgarcho <- function(y, p = 0.99){
  gfit <- garchFit(formula = ~garch(1, 1), data = y$oil,
                   cond.dist = "std", trace = FALSE)
  sigma <-  predict(gfit, n.ahead = 1)[3]
  df <- coef(gfit)["shape"]
  ES <- sigma * (dt(qt(p, df), df)/(1 - p)) *
    ((df + (qt(p, df))^2)/(df - 1))
  return(ES)
}
ESgarcho(goydatats,p=0.95)
##   standardDeviation
## 1          108.9479
ESgarchy <- function(y, p = 0.99){
  gfit <- garchFit(formula = ~garch(1, 1), data = y$yen,
                   cond.dist = "std", trace = FALSE)
  sigma <-  predict(gfit, n.ahead = 1)[3]
  df <- coef(gfit)["shape"]
  ES <- sigma * (dt(qt(p, df), df)/(1 - p)) *
    ((df + (qt(p, df))^2)/(df - 1))
  return(ES)
}
ESgarchy(goydatats,p=0.95)
## Warning in sqrt(diag(fit$cvar)): NaNs produced
##   standardDeviation
## 1           5.60357
goyreturns<- data.frame(mean(goydata$gold), mean(goydata$oil), mean(goydata$yen))
goyreturns
##   mean.goydata.gold. mean.goydata.oil. mean.goydata.yen.
## 1           550.2547          36.25545          160.3579
goycov<- cov(goyloss, use = "pairwise.complete.obs")
head(goycov)
##           gold       oil       yen
## gold 23.605351  6.809759 -3.399963
## oil   6.809759 70.396028  1.081998
## yen  -3.399963  1.081998  6.950542

Monthly continuously compounded returns as difference in log prices

ret.cc<-diff(log(goydatats))
portret<-apply(na.omit(ret.cc), 2, mean)
mu.hat.annual<-apply(na.omit(ret.cc), 2, mean) * 12   
cov.mat.annual<-cov(ret.cc,use = "pairwise.complete.obs") * 12
mu.hat.annual
##        gold         oil         yen 
##  0.07385320  0.05691500 -0.02445996
cov.mat.annual
##              gold         oil          yen
## gold  0.028326422 0.008171711 -0.004079956
## oil   0.008171711 0.084475234  0.001298398
## yen  -0.004079956 0.001298398  0.008340651
library(IntroCompFinR)
## Loading required package: xts
## Warning: package 'xts' was built under R version 3.4.4
## Loading required package: zoo
## Warning: package 'zoo' was built under R version 3.4.4
## 
## Attaching package: 'zoo'
## The following object is masked from 'package:timeSeries':
## 
##     time<-
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
GMVP<-globalMin.portfolio(mu.hat.annual,cov.mat.annual)
plot(GMVP)

GMVPweights<-GMVP$weights
wgold<-GMVPweights[1]
woil<-GMVPweights[2]
wyen<-GMVPweights[3]

Expected Shortfall measure of rik with 95% confidence.

lossg<-ESgarchg(goyloss,p=0.95)
losso<-ESgarcho(goyloss,p=0.95)
## Warning in sqrt(diag(fit$cvar)): NaNs produced
lossy<-ESgarchy(goyloss,p=0.95)
# Estimate GARCH model
# Step 1
gfitgoy<-lapply(goyloss,garchFit,formula=~garch(1,1),cond.dist="std",trace=FALSE)
## Warning in sqrt(diag(fit$cvar)): NaNs produced
gfitgoy
## $gold
## 
## Title:
##  GARCH Modelling 
## 
## Call:
##  FUN(formula = ..1, data = X[[i]], cond.dist = "std", trace = FALSE) 
## 
## Mean and Variance Equation:
##  data ~ garch(1, 1)
## <environment: 0x0000000020043610>
##  [data = X[[i]]]
## 
## Conditional Distribution:
##  std 
## 
## Coefficient(s):
##       mu     omega    alpha1     beta1     shape  
## -0.15603   0.61646   0.14811   0.83657   5.21538  
## 
## Std. Errors:
##  based on Hessian 
## 
## Error Analysis:
##         Estimate  Std. Error  t value Pr(>|t|)    
## mu      -0.15603     0.14687   -1.062 0.288086    
## omega    0.61646     0.37124    1.661 0.096809 .  
## alpha1   0.14811     0.04450    3.329 0.000873 ***
## beta1    0.83657     0.04469   18.719  < 2e-16 ***
## shape    5.21538     1.23216    4.233 2.31e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log Likelihood:
##  -1629.699    normalized:  -2.824434 
## 
## Description:
##  Wed Apr 03 01:02:49 2019 by user: Ravindra 
## 
## 
## $oil
## 
## Title:
##  GARCH Modelling 
## 
## Call:
##  FUN(formula = ..1, data = X[[i]], cond.dist = "std", trace = FALSE) 
## 
## Mean and Variance Equation:
##  data ~ garch(1, 1)
## <environment: 0x00000000266a8628>
##  [data = X[[i]]]
## 
## Conditional Distribution:
##  std 
## 
## Coefficient(s):
##          mu        omega       alpha1        beta1        shape  
## -1.4363e-05   7.0396e-05   1.0000e+00   3.7818e-01   2.8737e+00  
## 
## Std. Errors:
##  based on Hessian 
## 
## Error Analysis:
##          Estimate  Std. Error  t value Pr(>|t|)    
## mu     -1.436e-05   2.070e-03   -0.007    0.994    
## omega   7.040e-05          NA       NA       NA    
## alpha1  1.000e+00   4.786e-02   20.894   <2e-16 ***
## beta1   3.782e-01          NA       NA       NA    
## shape   2.874e+00   1.194e-01   24.076   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log Likelihood:
##  -1784.325    normalized:  -3.092418 
## 
## Description:
##  Wed Apr 03 01:02:49 2019 by user: Ravindra 
## 
## 
## $yen
## 
## Title:
##  GARCH Modelling 
## 
## Call:
##  FUN(formula = ..1, data = X[[i]], cond.dist = "std", trace = FALSE) 
## 
## Mean and Variance Equation:
##  data ~ garch(1, 1)
## <environment: 0x000000001df7e740>
##  [data = X[[i]]]
## 
## Conditional Distribution:
##  std 
## 
## Coefficient(s):
##      mu    omega   alpha1    beta1    shape  
## 0.10346  1.68901  0.14849  0.62233  6.71524  
## 
## Std. Errors:
##  based on Hessian 
## 
## Error Analysis:
##         Estimate  Std. Error  t value Pr(>|t|)    
## mu        0.1035      0.1072    0.965 0.334392    
## omega     1.6890      1.6576    1.019 0.308212    
## alpha1    0.1485      0.0777    1.911 0.055997 .  
## beta1     0.6223      0.2895    2.149 0.031603 *  
## shape     6.7152      2.0032    3.352 0.000802 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log Likelihood:
##  -1359.987    normalized:  -2.356997 
## 
## Description:
##  Wed Apr 03 01:02:50 2019 by user: Ravindra
mode(gfitgoy)
## [1] "list"
gprog<-unlist(lapply(gfitgoy,function(x) predict(x,n.ahead = 1)[3]))

# Estimate degrees-of-freedom parameters
gshape<-unlist(lapply(gfitgoy, function(x) x@fit$coef[5]))
# Can take a look at all paramaters of the GARCH model
gcoef<-unlist(lapply(gfitgoy, function(x) x@fit$coef))
# Estimates conditional standardized residuals are extracted.(h.t - conditional variance)
# Step 2
gresid<-as.matrix(data.frame(lapply(gfitgoy,function(x) x@residuals / sqrt(x@h.t))))
head(gresid)
##             gold          oil         yen
## [1,] -0.42024299 1.457163e-06  0.01103960
## [2,] -0.05372649 2.369522e-06 -0.03899564
## [3,] -0.04027726 3.853122e-06 -0.04238046
## [4,] -0.89690941 6.265604e-06 -0.03465949
## [5,]  0.27965462 1.018846e-05 -0.04677800
## [6,] -0.49908486 1.656689e-05 -0.04667386
#QQ plots of the standardized residuals
par(mfrow=c(2,2))
unlist(lapply(gfitgoy, function(x) plot(x, which=13)))
## $gold
## 
## Title:
##  GARCH Modelling 
## 
## Call:
##  FUN(formula = ..1, data = X[[i]], cond.dist = "std", trace = FALSE) 
## 
## Mean and Variance Equation:
##  data ~ garch(1, 1)
## <environment: 0x0000000020043610>
##  [data = X[[i]]]
## 
## Conditional Distribution:
##  std 
## 
## Coefficient(s):
##       mu     omega    alpha1     beta1     shape  
## -0.15603   0.61646   0.14811   0.83657   5.21538  
## 
## Std. Errors:
##  based on Hessian 
## 
## Error Analysis:
##         Estimate  Std. Error  t value Pr(>|t|)    
## mu      -0.15603     0.14687   -1.062 0.288086    
## omega    0.61646     0.37124    1.661 0.096809 .  
## alpha1   0.14811     0.04450    3.329 0.000873 ***
## beta1    0.83657     0.04469   18.719  < 2e-16 ***
## shape    5.21538     1.23216    4.233 2.31e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log Likelihood:
##  -1629.699    normalized:  -2.824434 
## 
## Description:
##  Wed Apr 03 01:02:49 2019 by user: Ravindra 
## 
## 
## $oil
## 
## Title:
##  GARCH Modelling 
## 
## Call:
##  FUN(formula = ..1, data = X[[i]], cond.dist = "std", trace = FALSE) 
## 
## Mean and Variance Equation:
##  data ~ garch(1, 1)
## <environment: 0x00000000266a8628>
##  [data = X[[i]]]
## 
## Conditional Distribution:
##  std 
## 
## Coefficient(s):
##          mu        omega       alpha1        beta1        shape  
## -1.4363e-05   7.0396e-05   1.0000e+00   3.7818e-01   2.8737e+00  
## 
## Std. Errors:
##  based on Hessian 
## 
## Error Analysis:
##          Estimate  Std. Error  t value Pr(>|t|)    
## mu     -1.436e-05   2.070e-03   -0.007    0.994    
## omega   7.040e-05          NA       NA       NA    
## alpha1  1.000e+00   4.786e-02   20.894   <2e-16 ***
## beta1   3.782e-01          NA       NA       NA    
## shape   2.874e+00   1.194e-01   24.076   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log Likelihood:
##  -1784.325    normalized:  -3.092418 
## 
## Description:
##  Wed Apr 03 01:02:49 2019 by user: Ravindra 
## 
## 
## $yen
## 
## Title:
##  GARCH Modelling 
## 
## Call:
##  FUN(formula = ..1, data = X[[i]], cond.dist = "std", trace = FALSE) 
## 
## Mean and Variance Equation:
##  data ~ garch(1, 1)
## <environment: 0x000000001df7e740>
##  [data = X[[i]]]
## 
## Conditional Distribution:
##  std 
## 
## Coefficient(s):
##      mu    omega   alpha1    beta1    shape  
## 0.10346  1.68901  0.14849  0.62233  6.71524  
## 
## Std. Errors:
##  based on Hessian 
## 
## Error Analysis:
##         Estimate  Std. Error  t value Pr(>|t|)    
## mu        0.1035      0.1072    0.965 0.334392    
## omega     1.6890      1.6576    1.019 0.308212    
## alpha1    0.1485      0.0777    1.911 0.055997 .  
## beta1     0.6223      0.2895    2.149 0.031603 *  
## shape     6.7152      2.0032    3.352 0.000802 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log Likelihood:
##  -1359.987    normalized:  -2.356997 
## 
## Description:
##  Wed Apr 03 01:02:50 2019 by user: Ravindra
#ACF of the squared residuals
par(mfrow=c(2,2))

unlist(lapply(gfitgoy, function(x) plot(x, which=11)))
## $gold
## 
## Title:
##  GARCH Modelling 
## 
## Call:
##  FUN(formula = ..1, data = X[[i]], cond.dist = "std", trace = FALSE) 
## 
## Mean and Variance Equation:
##  data ~ garch(1, 1)
## <environment: 0x0000000020043610>
##  [data = X[[i]]]
## 
## Conditional Distribution:
##  std 
## 
## Coefficient(s):
##       mu     omega    alpha1     beta1     shape  
## -0.15603   0.61646   0.14811   0.83657   5.21538  
## 
## Std. Errors:
##  based on Hessian 
## 
## Error Analysis:
##         Estimate  Std. Error  t value Pr(>|t|)    
## mu      -0.15603     0.14687   -1.062 0.288086    
## omega    0.61646     0.37124    1.661 0.096809 .  
## alpha1   0.14811     0.04450    3.329 0.000873 ***
## beta1    0.83657     0.04469   18.719  < 2e-16 ***
## shape    5.21538     1.23216    4.233 2.31e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log Likelihood:
##  -1629.699    normalized:  -2.824434 
## 
## Description:
##  Wed Apr 03 01:02:49 2019 by user: Ravindra 
## 
## 
## $oil
## 
## Title:
##  GARCH Modelling 
## 
## Call:
##  FUN(formula = ..1, data = X[[i]], cond.dist = "std", trace = FALSE) 
## 
## Mean and Variance Equation:
##  data ~ garch(1, 1)
## <environment: 0x00000000266a8628>
##  [data = X[[i]]]
## 
## Conditional Distribution:
##  std 
## 
## Coefficient(s):
##          mu        omega       alpha1        beta1        shape  
## -1.4363e-05   7.0396e-05   1.0000e+00   3.7818e-01   2.8737e+00  
## 
## Std. Errors:
##  based on Hessian 
## 
## Error Analysis:
##          Estimate  Std. Error  t value Pr(>|t|)    
## mu     -1.436e-05   2.070e-03   -0.007    0.994    
## omega   7.040e-05          NA       NA       NA    
## alpha1  1.000e+00   4.786e-02   20.894   <2e-16 ***
## beta1   3.782e-01          NA       NA       NA    
## shape   2.874e+00   1.194e-01   24.076   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log Likelihood:
##  -1784.325    normalized:  -3.092418 
## 
## Description:
##  Wed Apr 03 01:02:49 2019 by user: Ravindra 
## 
## 
## $yen
## 
## Title:
##  GARCH Modelling 
## 
## Call:
##  FUN(formula = ..1, data = X[[i]], cond.dist = "std", trace = FALSE) 
## 
## Mean and Variance Equation:
##  data ~ garch(1, 1)
## <environment: 0x000000001df7e740>
##  [data = X[[i]]]
## 
## Conditional Distribution:
##  std 
## 
## Coefficient(s):
##      mu    omega   alpha1    beta1    shape  
## 0.10346  1.68901  0.14849  0.62233  6.71524  
## 
## Std. Errors:
##  based on Hessian 
## 
## Error Analysis:
##         Estimate  Std. Error  t value Pr(>|t|)    
## mu        0.1035      0.1072    0.965 0.334392    
## omega     1.6890      1.6576    1.019 0.308212    
## alpha1    0.1485      0.0777    1.911 0.055997 .  
## beta1     0.6223      0.2895    2.149 0.031603 *  
## shape     6.7152      2.0032    3.352 0.000802 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log Likelihood:
##  -1359.987    normalized:  -2.356997 
## 
## Description:
##  Wed Apr 03 01:02:50 2019 by user: Ravindra

## Copula
# pseudo-uniform variables that generates probabilites for each risk 
# (measured as conditional resid)
# Step 3
U <- sapply(1:3, function(y) pt(gresid[, y], df = gshape[y]))
head(U)
##           [,1]      [,2]      [,3]
## [1,] 0.3455316 0.5000005 0.5042437
## [2,] 0.4795754 0.5000009 0.4850139
## [3,] 0.4846844 0.5000014 0.4837139
## [4,] 0.2046283 0.5000023 0.4866794
## [5,] 0.6047519 0.5000037 0.4820254
## [6,] 0.3190293 0.5000061 0.4820654
hist(U)

# Student's t copula is estimated based on Kendall's rank correlations.  
# Step 4
cop <- fit.tcopula(Udata = U, method = "Kendall")
## Warning in FUN(newX[, i], ...): NaNs produced

## Warning in FUN(newX[, i], ...): NaNs produced

## Warning in FUN(newX[, i], ...): NaNs produced
## Warning in dmt(Qdata, df, rep(0, d), Sigma, log = TRUE): value out of range
## in 'lgamma'

## Warning in dmt(Qdata, df, rep(0, d), Sigma, log = TRUE): value out of range
## in 'lgamma'
## Warning in log(pi * df): NaNs produced
## Warning in nlminb(startdf, negloglik2, data = Udata, P = P, ...): NA/NaN
## function evaluation
## Warning in FUN(newX[, i], ...): NaNs produced

## Warning in FUN(newX[, i], ...): NaNs produced

## Warning in FUN(newX[, i], ...): NaNs produced
## Warning in dmt(Qdata, df, rep(0, d), Sigma, log = TRUE): value out of range
## in 'lgamma'

## Warning in dmt(Qdata, df, rep(0, d), Sigma, log = TRUE): value out of range
## in 'lgamma'
## Warning in log(pi * df): NaNs produced
## Warning in nlminb(startdf, negloglik2, data = Udata, P = P, ...): NA/NaN
## function evaluation
# 100,000 random losses simulated for each financial instrument 
# Step 5
rcop <- rcopula.t(100000, df = cop$nu, Sigma = cop$P)
head(rcop)
##           [,1]        [,2]       [,3]
## [1,] 0.7410856 0.904027185 0.10578590
## [2,] 0.1736152 0.172181024 0.19423294
## [3,] 0.5297456 0.803959923 0.38045824
## [4,] 0.4972565 0.183309712 0.58877831
## [5,] 0.2278439 0.850448726 0.11266660
## [6,] 0.9942562 0.003801331 0.02126613
hist(rcop)

# Compute the quantiles for these Monte Carlo draws.
#Step 6
qcop <- sapply(1:3, function(x) qstd(rcop[, x], nu = gshape[x]))
head(qcop)
##              [,1]       [,2]       [,3]
## [1,]  0.544306887  0.9366054 -1.1569664
## [2,] -0.810869788 -0.6223549 -0.7725735
## [3,]  0.061476644  0.5538579 -0.2655790
## [4,] -0.005663318 -0.5893127  0.1955359
## [5,] -0.632445123  0.6968221 -1.1187298
## [6,]  2.995385310 -3.7212698 -2.0928896
hist(qcop)

# creating a matix of 1 period ahead predictions of standard deviations
ht.mat <- matrix(gprog, nrow = 100000, ncol = ncol(goyloss), byrow = TRUE)
head(ht.mat)
##          [,1]     [,2]     [,3]
## [1,] 2.900298 10.19218 2.375486
## [2,] 2.900298 10.19218 2.375486
## [3,] 2.900298 10.19218 2.375486
## [4,] 2.900298 10.19218 2.375486
## [5,] 2.900298 10.19218 2.375486
## [6,] 2.900298 10.19218 2.375486
pf <- qcop * ht.mat
head(pf)
##             [,1]       [,2]       [,3]
## [1,]  1.57865239   9.546050 -2.7483569
## [2,] -2.35176434  -6.343153 -1.8352372
## [3,]  0.17830061   5.645019 -0.6308791
## [4,] -0.01642531  -6.006381  0.4644928
## [5,] -1.83427957   7.102136 -2.6575264
## [6,]  8.68751120 -37.927848 -4.9716288
## ES 95 percent
weights <- c(wgold,woil,wyen)
# The simulated portfolio losses are then determined 
# as the outcome of the matrix-weight vector product
# Step 7
pfall <- (qcop * ht.mat) %*% weights
head(pfall)
##            [,1]
## [1,] -1.3334048
## [2,] -2.0671034
## [3,] -0.2839216
## [4,]  0.2015262
## [5,] -2.2346537
## [6,] -1.9911827
# Step 8
pfall.es95 <- median(tail(sort(pfall), 5000))
pfall.es95
## [1] 3.316646
pfall.var95 <- min(tail(sort(pfall), 5000))
pfall.var95
## [1] 2.480332