CoCalc -- Risk_budget

GitHub Repository: braverock/portfolioanalytics
Path: blob/master/sandbox/riskbudgetpaper(superseded)/R_Allocation/Risk_budget_functions.R
¹⁴³³ views
1

2
PortfolioOptim = function( minriskcriterion = "mES" , MinMaxComp = F, percriskcontribcriterion = "mES" , 
3
                     R = NULL, mu = NULL , sigma = NULL, M3=NULL,M4=NULL, alpha = 0.05, alphariskbudget = 0.05,
4
                     lower = NULL , upper = NULL, Riskupper = NULL, Returnlower = NULL, RBlower = NULL , RBupper = NULL, precision = 1e-3 , 
5
                     controlDE = NULL , optimize_method = "DEoptim+L-BFGS-B", penalty = NULL ){
6

7
     # Description: 
8
     # This function produces the portfolio that minimimizes 
9
     # either the portfolio risk (when MinMaxComp = F) or the portfolio concentration (when MinMaxComp = T) 
10
     # subject to
11
     # lower <= weights <= upper
12
     # risk <= Riskupper
13
     # expected return >= Returnlower 
14
     # RBlower <= percentage risk <= RBupper
15

16
     # Input: 
17
     # Either the multivariate return series is given or estimates of the mean, covariance, coskewness or cokurtosis
18
     require(zoo); require(PerformanceAnalytics);   require(DEoptim)
19

20
     if( !is.null(R) ){
21
          R = clean.boudt2( R , alpha = alphariskbudget )[[1]];
22
          T = nrow(R);N=ncol(R)
23
          mu = matrix( as.vector(apply(R,2,'mean')),ncol=1);
24
          sigma = cov(R);
25
          M3 = PerformanceAnalytics:::M3.MM(R)
26
          M4 = PerformanceAnalytics:::M4.MM(R)
27
     }else{ N = length(mu) }
28

29
     if( is.null(lower)     ){ lower = rep(0,N)      }  ; if( is.null(upper)       ){ upper = rep(1,N)     }
30
     if( is.null(RBlower)   ){ RBlower = rep(-Inf,N) }  ; if( is.null(RBupper)     ){ RBupper = rep(Inf,N) }
31
     if( is.null(Riskupper) ){ Riskupper = Inf       }  ; if( is.null(Returnlower) ){ Returnlower = -Inf   }  
32

33
     switch( percriskcontribcriterion ,
34
                StdDev = { percriskcontrib = function(w){ return( Portsd(w,mu=mu,sigma=sigma, precision=9)                                           ) }},
35
                GVaR   = { percriskcontrib = function(w){ return( PortgausVaR(w,alpha=alphariskbudget,mu=mu,sigma=sigma, precision=9)                ) }},
36
                GES    = { percriskcontrib = function(w){ return( PortgausES(w,mu=mu,alpha=alphariskbudget,sigma=sigma, precision=9)                 ) }},
37
                mVaR   = { percriskcontrib = function(w){ return( PortMVaR(w,mu=mu,alpha=alphariskbudget,sigma=sigma,M3=M3,M4=M4, precision=9)       ) }},
38
                mES    = { percriskcontrib = function(w){ return( operPortMES(w,mu=mu,alpha=alphariskbudget,sigma=sigma,M3=M3,M4=M4, precision=9)    ) }}
39
     ) #end function that finds out which percentage risk contribution criterion to use
40
     switch( minriskcriterion ,
41
                StdDev = { prisk = function(w){ return( stddevfun(w,mu=mu,sigma=sigma) ) }},
42
                GVaR   = { prisk = function(w){ return( gausVaRfun(w,alpha=alpha,mu=mu,sigma=sigma) ) }},
43
                GES    = { prisk = function(w){ return( gausESfun(w,mu=mu,alpha=alpha,sigma=sigma) )  }},
44
                mVaR   = { prisk = function(w){ return( MVaRfun(w,mu=mu,alpha=alpha,sigma=sigma,M3=M3,M4=M4) )}},
45
                mES    = { prisk = function(w){ return( operMESfun(w,mu=mu,alpha=alpha,sigma=sigma,M3=M3,M4=M4)) }}
46
     ) #end function that finds out which risk function to minimize
47

48
    if( is.null(penalty) ){ 
49
       abspenalty = 1e6; relpenalty = 100*N 
50
    }else{
51
       abspenalty = relpenalty = penalty; 
52
    }; 
53
    if( optimize_method == "DEoptim" | optimize_method == "DEoptim+L-BFGS-B" ){
54
        objective = function( w ){
55
           w = w/sum(w)  
56
           cont = percriskcontrib( w );  percrisk = cont[[3]]; crisk = cont[[2]] ; 
57
           if(MinMaxComp){ out = max( crisk ) }else{ out = prisk(w) }
58
           # add weight constraints
59
           out_con = out + out*abspenalty*sum( 1*( w < lower )+1*( w > upper ) )
60
           # add risk budget constraint 
61
           con1 = sum(  (percrisk-RBupper)*( percrisk > RBupper  ),na.rm=TRUE ) ; if( is.na(con1) ){ con1 = 0 } # because of Inf*0
62
           con2 = sum(  (RBlower-percrisk)*( percrisk < RBlower  ),na.rm=TRUE  ); if( is.na(con2) ){ con2 = 0 }  
63
           out_con = out_con + out*relpenalty*(con1+con2)
64
           # add minimum risk constraint
65
           con = ( prisk(w) - Riskupper)*( prisk(w) > Riskupper ); if( is.na(con) ){ con = 0 }  
66
           out_con = out_con + out*relpenalty*con
67
           # add minimum return constraint 
68
           # portfolio return and risk are in the same unit, but portfolio return is typically some orders of magnitude smaller
69
           # say: as a very conservative choice: 100 
70
           preturn = sum( w*mu ) ;
71
           con = ( Returnlower - preturn)*( preturn < Returnlower ); if( is.na(con) ){ con = 0 }
72
           out_con = out_con + out*100*relpenalty*con
73
           return(out_con)
74
         }
75

76
         if( is.null(controlDE) ){
77
            # initial population is generated with random_portfolios function
78
            require('PortfolioAnalytics')
79
            eps = 0.01
80
            rpconstraint<-constraint(assets=length(lower), min_sum=(1-eps), max_sum=(1+eps), 
81
              min=lower, max=upper, weight_seq=generatesequence())
82
            rp<- random_portfolios(rpconstraints=rpconstraint,permutations=200)
83
            controlDE = list(  NP=as.numeric(nrow(rp)) , initialpop=rp,trace=F ) 
84
         }
85
         # print(controlDE)
86
         minw = DEoptim( fn = objective ,  lower = lower , upper = upper , control = controlDE)
87
         fvalues = minw$member$bestval
88
         diff = as.vector( quantile(fvalues,0.10) - min(fvalues) )
89
         print(c("diff",diff))
90
         best = min( fvalues ) ; print(best)
91
         bestsol = minw ; 
92

93
         while( diff>1e-6 ){
94
            pop = as.matrix(minw$member$pop)
95
            pop[1,] = minw$optim$bestmem;
96
            minw = DEoptim( fn = objective ,  lower = lower , upper = upper , 
97
                control = list( itermax = 150, NP=as.numeric(nrow(pop)) , initialpop=pop,trace=F ) )
98
            fvalues = minw$member$bestval
99
            diff = best -  min(fvalues)  ; 
100
            if( diff > 0 ){ best = min( fvalues ) ; bestsol = minw ; print(best) }
101
          }
102
          minw = bestsol$optim$bestmem/sum(bestsol$optim$bestmem)  ; #full investment constraint
103
    }
104

105

106
    if( optimize_method == "DEoptim+L-BFGS-B" ){
107
         print(minw) ; print( objective(minw) )
108
         print("local improvement")
109
         out = optim( par=minw, f = objective, lower = lower, upper = upper, method="L-BFGS-B" ) 
110
         minw = out$par/sum(out$par)
111
         print(minw) ; print( objective(minw) )
112
    }
113

114
    if( optimize_method == "L-BFGS-B" ){
115
        objective = function( w ){
116
           if(sum(w)==0){w=w+0.001}
117
           w = w/sum(w)  
118
           cont = percriskcontrib( w );  percrisk = cont[[3]]; crisk = cont[[2]] ; 
119
           if(MinMaxComp){ out = max( crisk ) }else{ out = prisk(w) }
120
           # add risk budget constraint 
121
           con1 = sum(  (percrisk-RBupper)*( percrisk > RBupper  ),na.rm=TRUE ) ; if( is.na(con1) ){ con1 = 0 } # because of Inf*0
122
           con2 = sum(  (RBlower-percrisk)*( percrisk < RBlower  ),na.rm=TRUE  ); if( is.na(con2) ){ con2 = 0 }  
123
           out_con = out + out*relpenalty*(con1+con2)
124
           # add minimum risk constraint
125
           con = ( prisk(w) - Riskupper)*( prisk(w) > Riskupper ); if( is.na(con) ){ con = 0 }  
126
           out_con = out_con + out*relpenalty*con
127
           return(out_con)
128
         }
129

130
         if(Returnlower==-Inf){ 
131
            inittheta = rep(1,N)/N 
132
            out = optim( par=inittheta, f = objective, lower = lower, upper = upper ) 
133
         }else{
134
            Amat = rbind(diag(x =1,nrow=N,ncol=N), diag(x =-1,nrow=N,ncol=N), rep(1,N), rep(-1,N),as.vector(mu))
135
            inittheta = rep(0.001,N);
136
            inittheta[mu==max(mu)] = 1; inittheta = 1-sum(inittheta[mu!=max(mu)]  );   
137
            out = constrOptim( theta=inittheta, f = objective, grad=NULL,ui=Amat,
138
                            ci = c(rep(0,N),rep(-1,N),0.99999,-1.0001,Returnlower) ) 
139
         }
140

141
         minw = out$par/sum(out$par)
142
    }
143
    cont = percriskcontrib( minw );  percrisk = cont[[3]];    crisk = cont[[2]] ;    
144
    # check 
145
    print( "out = list( minw , sum( minw*mu ) , prisk(minw) , percriskcontrib(minw)" )
146
    out = list( minw , sum( minw*mu ) , prisk(minw) , percrisk , crisk )
147
    print( out )
148
    return(out)
149
}
150

151
findportfolio.dynamic = function(R, from, to, names.input = NA,  names.assets, 
152
              	p = 0.95 , priskbudget = 0.95,  mincriterion = "mES" , percriskcontribcriterion = "mES"  ,
153
                  strategy , controlDE =  list( VTR = 0 , NP=200 , itermax = 200,trace=F ), optimize_method = "DEoptim+L-BFGS-B" )
154
{ # @author Kris Boudt and Brian G. Peterson
155

156
    # Description:
157
    #
158
    # Performs a loop over the reallocation periods with estimation samples given by from:to
159
    # It calls the function RBconportfolio to obtain the optimal weights of the strategy.
160
    #
161
    # @todo
162
    #
163
    # R                 matrix/zoo holding historical returns on risky assets
164
    #
165
    # names             vector holding the names of the .csv files to be read
166
    #
167
    # from, to          define the estimation sample
168
    #
169
    # criteria	      the criterion to be optimized
170
    #
171
    # columns.crit      the columns of R in which the criteria are located
172
    #
173
    # percriskcontribcriterion risk measure used for the risk budget constraints
174
    #
175
    # strategy = c( "EqualRisk" , "EqualWeight" , "MinRisk" , "MinRiskConc" , 
176
    #                "MinRisk_PositionLimit"     , "MinRisk_RiskLimit"     , "MinRisk_ReturnTarget",
177
    #                "MinRiskConc_PositionLimit" , "MinRiskConc_RiskLimit" , "MinRiskConc_ReturnTarget")
178

179
    # Return:
180
    # List with first element optimal weights per reallocation period and associated percentage CVaR contributions. 
181

182
    # Create a matrix that will hold for each method and each vector the best weights
183

184
    cPeriods = length(from);    
185

186
    out = matrix(  rep(0, cPeriods*(cAssets)) , ncol= (cAssets) );
187
    RCout = matrix(  rep(0, cPeriods*(cAssets)) , ncol= (cAssets) );
188
    RiskReturnout = matrix(  rep(0, cPeriods*2 ) , ncol= 2 );
189
 
190
    # first cPeriods rows correspond to cCriteria[1] and so on
191

192
    # downside risk
193
    alpha = 1 - p;
194
    alphariskbudget = 1-priskbudget;
195

196
    # Estimation of the return mean vector, covariance, coskewness and cokurtosis matrix
197

198

199
    if(strategy=="EqualRisk"){
200
       lower = rep(0,cAssets); upper=rep(1,cAssets) 
201
       RBlower = rep(1/cAssets,cAssets) ; RBupper = rep(1/cAssets,cAssets)  ;
202
    }
203

204
    if(strategy=="EqualWeight"){
205
       lower = rep(1/cAssets,cAssets); upper=rep(1/cAssets,cAssets) 
206
       RBlower = rep(-Inf,cAssets) ; RBupper = rep(Inf,cAssets)  ;
207
    }
208

209
    if(strategy=="MinRisk" | strategy=="MinRiskConc" | strategy=="MinRisk_ReturnTarget" | strategy=="MinRiskConc_ReturnTarget"){
210
       lower = rep(0,cAssets); upper=rep(1,cAssets) 
211
       RBlower = rep(-Inf,cAssets) ; RBupper = rep(Inf,cAssets)  ;
212
    }
213

214
    MinMaxComp = F; mutarget = -Inf; 
215
    if( strategy=="MinRiskConc" | strategy=="MinRiskConc_PositionLimit" | strategy=="MinRiskConc_RiskLimit" | strategy=="MinRiskConc_ReturnTarget"  ){
216
        MinMaxComp = T;
217
    }
218

219
   if(strategy=="MinRisk_PositionLimit" | strategy=="MinRiskConc_PositionLimit"){
220
      lower = rep(0,cAssets); upper=rep(0.4,cAssets) 
221
      RBlower = rep(-Inf,cAssets) ; RBupper = rep(Inf,cAssets)  ;
222
   }
223

224
   if(strategy=="MinRisk_RiskLimit" | strategy=="MinRiskConc_RiskLimit"){
225
      lower = rep(0,cAssets); upper=rep(1,cAssets) 
226
      RBlower = rep(-Inf,cAssets) ; RBupper = rep(0.40,cAssets)  ;
227
   }
228

229
   # initial population is generated with random_portfolios function
230
   require('PortfolioAnalytics')
231
   eps = 0.01
232
   rpconstraint<-constraint(assets=length(lower), min_sum=(1-eps), max_sum=(1+eps), 
233
             min=lower, max=upper, weight_seq=generatesequence())
234
   rp<- random_portfolios(rpconstraints=rpconstraint,permutations=controlDE$NP)
235
   controlDE = list(  NP=as.numeric(nrow(rp)) , initialpop=rp,trace=F ) 
236

237
    for( per in c(1:cPeriods) ){
238

239
       print("-----------New estimation period ends on observation------------------")
240
       print( paste(to[per],"out of total number of obs equal to", max(to) ));
241
       print("----------------------------------------------------------------")
242

243
       # Estimate GARCH model with data from inception
244

245
       inception.R = window(R, start = as.Date(from[1]) , end = as.Date(to[per]) );
246

247
       # Estimate comoments of innovations with rolling estimation windows
248
       in.sample.R = window(R, start = as.Date(from[per]) , end = as.Date(to[per]) );
249
       in.sample.R = checkData(in.sample.R, method="matrix"); 
250

251
       # Estimation of mean return
252
       M = c();
253
       library(TTR)
254
       Tmean = 47 # monthly returns: 4 year exponentially weighted moving average
255
       for( i in 1:cAssets ){
256
         M = cbind( M , as.vector( EMA(x=inception.R[,i],n=Tmean) ) ) #2/(n+1)
257
       }
258
       M = zoo( M , order.by=time(inception.R) )
259

260
       # Center returns (shift by one observations since M[t,] is rolling mean t-Tmean+1,...,t; otherwise lookahead bias)
261
       inception.R.cent = inception.R;
262
       ZZ = matrix( rep(as.vector( apply( inception.R[1:Tmean, ] , 2 , 'mean' )),Tmean),byrow=T,nrow=Tmean);
263
       inception.R.cent[1:Tmean,] = inception.R[1:Tmean, ] - ZZ;
264
       if( nrow(inception.R)>(Tmean+1) ){ 
265
                 A = M[Tmean:(nrow(inception.R)-1),];
266
                 A = zoo( A , order.by = time(inception.R[(Tmean+1):nrow(inception.R), ])) ; #shift dates; otherwise zoo poses problem
267
                 inception.R.cent[(Tmean+1):nrow(inception.R), ] = inception.R[(Tmean+1):nrow(inception.R), ] - A}
268
       
269
       # Garch estimation 
270
       S = c();
271
       for( i in 1:cAssets ){
272
            gout =  garchFit(formula ~ garch(1,1), data = inception.R.cent[,i],include.mean = F, cond.dist="QMLE", trace = FALSE )
273
            if( as.vector(gout@fit$coef["alpha1"]) < 0.01 ){
274
               sigmat = rep( sd( as.vector(inception.R.cent[,i])), length(inception.R.cent[,i]) ); 
275
            }else{
276
               sigmat = gout@sigma.t
277
            }
278
            S = cbind( S , sigmat)
279
       }
280
       S = zoo( S , order.by=time(inception.R.cent) )
281

282
       # Estimate correlation, coskewness and cokurtosis matrix locally using cleaned innovation series in three year estimation window
283
       selectU = window(inception.R.cent, start = as.Date(from[per]) , end = as.Date(to[per]) )
284
       selectU = selectU/window(S, start = as.Date(from[per]) , end = as.Date(to[per]) );
285
       selectU = clean.boudt2(selectU , alpha = 0.05 )[[1]];
286
       Rcor = cor(selectU)
287
       D = diag( as.vector(tail(S,n=1)  ),ncol=cAssets )
288
       sigma = D%*%Rcor%*%D
289

290
       # we only need mean and conditional covariance matrix of last observation
291
       mu = matrix(tail(M,n=1),ncol=1 ) ;
292
       D = diag( as.vector(as.vector(tail(S,n=1) ) ),ncol=cAssets )
293
       sigma = D%*%Rcor%*%D
294
       in.sample.T = nrow(selectU);
295
       # set volatility of all U to last observation, such that cov(rescaled U)=sigma 
296
       selectU = selectU*matrix( rep(as.vector(tail(S,n=1)),in.sample.T  ) , ncol = cAssets , byrow = T )
297
       M3 = PerformanceAnalytics:::M3.MM(selectU)
298
       M4 = PerformanceAnalytics:::M4.MM(selectU)
299

300
       
301
       mESfun = function(series){ return( operMES(series,alpha=alpha,2) ) }
302
       
303
 
304
       if(strategy=="MinRisk_ReturnTarget" | strategy=="MinRiskConc_ReturnTarget"){
305
           mutarget = mean( mu );
306
           print( c("Minimum return requirement is" , mutarget) )
307
       }
308

309
       if(strategy=="EqualWeight"){
310
           sol1 = rep(1/cAssets,cAssets);
311
           switch( percriskcontribcriterion ,
312
                StdDev = { percriskcontrib = function(w){ cont = Portsd(w,mu=mu,sigma=sigma) ; return( cont  ) }},
313
                GVaR = {percriskcontrib = function(w){ cont =  PortgausVaR(w,alpha=alphariskbudget,mu=mu,sigma=sigma) ; return( cont ) }},
314
                GES = {percriskcontrib = function(w){ cont =  PortgausES(w,mu=mu,alpha=alphariskbudget,sigma=sigma) ; return( cont  ) }},
315
                mVaR = {percriskcontrib = function(w){ cont = PortMVaR(w,mu=mu,alpha=alphariskbudget,sigma=sigma,M3=M3,M4=M4) ; return( cont  ) }},
316
                mES = {percriskcontrib = function(w){ cont = operPortMES(w,mu=mu,alpha=alphariskbudget,sigma=sigma,M3=M3,M4=M4) ; return( cont ) }
317
                }
318
           )
319
           sol2 = percriskcontrib( sol1 )
320
           solution = list( sol1 , sol2[[3]] , sol2[[1]] , mean(mu) ) ;   
321
       }else{
322
           solution = PortfolioOptim(  minriskcriterion = "mES" , MinMaxComp = MinMaxComp, percriskcontribcriterion = "mES" , 
323
                                       mu = mu , sigma = sigma, M3=M3 , M4=M4 , alpha = alpha , alphariskbudget = alphariskbudget , 
324
                                       lower = lower , upper = upper , Riskupper = Inf , Returnlower= mutarget , RBlower = RBlower, RBupper = RBupper , 
325
                                       controlDE = controlDE , optimize_method = optimize_method ) 
326
           # [[1]] : weights ; [[2]] : expected portfolio return ; [[3]] : portfolio CVaR
327
           # [[4]] : percentage CVaR ; [[5]] : component CVaR
328
           solution = list( solution[[1]] , solution[[4]] , solution[[3]] , solution[[2]] ); 
329
       }
330
       out[ per, ]    = as.vector( solution[[1]] ) 
331
       RCout[per,  ]  = as.vector( solution[[2]] )   
332
       RiskReturnout[per,] = c( as.numeric(solution[[3]]) , as.numeric(solution[[4]]) )
333

334
 }#end loop over the rebalancing periods; indexed by per=1,...,cPeriods
335

336
  # Output save
337
  rownames(out) = rownames(RCout) = rownames(RiskReturnout) =  names.input; 
338
  colnames(out) = colnames(RCout) = names.assets;
339
  colnames(RiskReturnout) = c( "risk" , "expreturn" )
340

341
  EWweights = c( rep(1/cAssets,cAssets) )
342
  EWweights = matrix ( rep(EWweights,cPeriods) , ncol=(cAssets) , byrow = TRUE )
343
  rownames(EWweights) = names.input; colnames(EWweights) = names.assets;
344

345
  return( list( out, RCout, RiskReturnout) ) 
346
}
347

348
clean.boudt2 = 
349
function (R, alpha = 0.01, trim = 0.001) 
350
{
351
    # @author Kris Boudt and Brian Peterson
352
    # Cleaning method as described in 
353
    #     Boudt, Peterson and Croux. 2009. Estimation and decomposition of downside risk for portfolios with non-normal returns. Journal of Risk.
354

355
    stopifnot("package:robustbase" %in% search() || require("robustbase", 
356
        quietly = TRUE))
357
    R = checkData(R, method = "zoo")
358
    T = dim(R)[1]
359
    date = c(1:T)
360
    N = dim(R)[2]
361
    MCD = covMcd(as.matrix(R), alpha = 1 - alpha)
362
    # mu = as.matrix(MCD$raw.center)
363
    mu = MCD$raw.center
364
    sigma = MCD$raw.cov
365
    invSigma = solve(sigma)
366
    vd2t = c()
367
    cleaneddata = R
368
    outlierdate = c()
369
    for (t in c(1:T)) {
370
        d2t = as.matrix(R[t, ] - mu) %*% invSigma %*% t(as.matrix(R[t,] - mu))
371
        vd2t = c(vd2t, d2t)
372
    }
373
    out = sort(vd2t, index.return = TRUE)
374
    sortvd2t = out$x
375
    sortt = out$ix
376
    empirical.threshold = sortvd2t[floor((1 - alpha) * T)]
377
    T.alpha = floor(T * (1 - alpha)) + 1
378
    cleanedt = sortt[c(T.alpha:T)]
379
    for (t in cleanedt) {
380
        if (vd2t[t] > qchisq(1 - trim, N)) {
381
            cleaneddata[t, ] = sqrt(max(empirical.threshold, 
382
                qchisq(1 - trim, N))/vd2t[t]) * R[t, ]
383
            outlierdate = c(outlierdate, date[t])
384
        }
385
    }
386
    return(list(cleaneddata, outlierdate))
387
}
388

389

390
Ipower = function(power,h)
391
{
392

393
   fullprod = 1;
394

395
   if( (power%%2)==0 ) #even number: number mod is zero
396
   {
397
      pstar = power/2;
398
      for(j in c(1:pstar))
399
      {
400
         fullprod = fullprod*(2*j)
401
      }
402
      I = fullprod*dnorm(h);
403

404
      for(i in c(1:pstar) )
405
      { 
406
         prod = 1;
407
         for(j in c(1:i) )
408
         {
409
            prod = prod*(2*j)
410
         }
411
         I = I + (fullprod/prod)*(h^(2*i))*dnorm(h)
412
      }
413
   }
414
   else{
415
      pstar = (power-1)/2
416
      for(j in c(0:pstar) )
417
      {
418
         fullprod = fullprod*( (2*j)+1 )
419
      }
420
      I = -fullprod*pnorm(h);
421

422
      for(i in c(0:pstar) )
423
      { 
424
         prod = 1;
425
         for(j in c(0:i) )
426
         {
427
            prod = prod*( (2*j) + 1 )
428
         }
429
         I = I + (fullprod/prod)*(h^(  (2*i) + 1))*dnorm(h)
430
      }
431
   }
432
   return(I)
433
}
434

435

436
# Definition of statistics needed to compute Gaussian and modified VaR and ES for the return series of portfolios
437
# and to compute the contributions to portfolio downside risk, made by the different positions in the portfolio.
438
#----------------------------------------------------------------------------------------------------------------
439

440

441

442
m2 = function(w,sigma)
443
{
444
   return(t(w)%*%sigma%*%w)
445
}
446
derm2 = function(w,sigma)
447
{
448
   return(2*sigma%*%w)
449
}
450
m3 = function(w,M3)
451
{
452
   return(t(w)%*%M3%*%(w%x%w))
453
}
454
derm3 = function(w,M3)
455
{
456
   return(3*M3%*%(w%x%w))
457
}
458
m4 = function(w,M4)
459
{
460
   return(t(w)%*%M4%*%(w%x%w%x%w))
461
}
462
derm4 = function(w,M4)
463
{
464
   return(4*M4%*%(w%x%w%x%w))
465
}
466

467
StdDevfun = function(w,sigma){ return(  sqrt( t(w)%*%sigma%*%w  )) }
468

469
GVaRfun = function(w,alpha,mu,sigma){ return (- (t(w)%*%mu) - qnorm(alpha)*sqrt( t(w)%*%sigma%*%w  ) ) }
470

471
mVaRfun = function(w,alpha,mu,sigma,M3,M4){
472
    pm4 = t(w)%*%M4%*%(w%x%w%x%w) ; pm3 = t(w)%*%M3%*%(w%x%w) ; pm2 =  t(w)%*%sigma%*%w ; 
473
    skew = pm3 / pm2^(3/2);
474
    exkurt = pm4 / pm2^(2) - 3; z = qnorm(alpha);
475
    h = z + (1/6)*(z^2 -1)*skew
476
    h = h + (1/24)*(z^3 - 3*z)*exkurt - (1/36)*(2*z^3 - 5*z)*skew^2;
477
    return (- (t(w)%*%mu) - h*sqrt( pm2  ) ) }
478

479
resmVaRfun = function(w,alpha,mu,sigma,ressigma,M3,M4){
480
    pm4 = t(w)%*%M4%*%(w%x%w%x%w) ; pm3 = t(w)%*%M3%*%(w%x%w) ; pm2 =  t(w)%*%sigma%*%w ; respm2 =  t(w)%*%resSigma%*%w ;
481
    skew = pm3 / respm2^(3/2);
482
    exkurt = pm4 / respm2^(2) - 3; z = qnorm(alpha);
483
    h = z + (1/6)*(z^2 -1)*skew
484
    h = h + (1/24)*(z^3 - 3*z)*exkurt - (1/36)*(2*z^3 - 5*z)*skew^2;
485
    return (- (t(w)%*%mu) - h*sqrt( pm2  ) ) }
486

487
GESfun = function(w,alpha,mu,sigma,M3,M4){
488
    return (- (t(w)%*%mu) + dnorm(qnorm(alpha))*sqrt(t(w)%*%sigma%*%w)/alpha ) }
489

490
operMESfun = function(w,alpha,mu,sigma,M3,M4){
491
    pm4 = t(w)%*%M4%*%(w%x%w%x%w) ; pm3 = t(w)%*%M3%*%(w%x%w) ; pm2 =  t(w)%*%sigma%*%w ; 
492
    skew = pm3 / pm2^(3/2);
493
    exkurt = pm4 / pm2^(2) - 3; z = qnorm(alpha);
494
    h = z + (1/6)*(z^2 -1)*skew
495
    h = h + (1/24)*(z^3 - 3*z)*exkurt - (1/36)*(2*z^3 - 5*z)*skew^2;
496
    E = dnorm(h)
497
    E = E + (1/24)*(   Ipower(4,h) - 6*Ipower(2,h) + 3*dnorm(h)   )*exkurt
498
    E = E +  (1/6)*(   Ipower(3,h) - 3*Ipower(1,h)   )*skew;
499
    E = E + (1/72)*(  Ipower(6,h) -15*Ipower(4,h)+ 45*Ipower(2,h) - 15*dnorm(h) )*(skew^2)
500
    E = E/alpha
501
    return (- (t(w)%*%mu) - sqrt(pm2)*min(-E,h) ) }
502

503

504
Portmean = function(w,mu,precision=4)
505
{
506
   return( list(  round( t(w)%*%mu , precision) , round ( as.vector(w)*as.vector(mu) , precision ) , round( as.vector(w)*as.vector(mu)/t(w)%*%mu) , precision) )
507
}
508

509
Portsd =  function(w,sigma,precision=4)
510
{
511
   pm2 = m2(w,sigma)
512
   dpm2 = derm2(w,sigma)
513
   dersd = (0.5*as.vector(dpm2))/sqrt(pm2);
514
   contrib = dersd*as.vector(w)
515
   return(list(  round( sqrt(pm2) , precision ) , round( contrib , precision ) , round ( contrib/sqrt(pm2) , precision) )) 
516
}
517

518

519
PortgausVaR =  function(alpha,w,mu,sigma,precision=4){
520
   location = t(w)%*%mu
521
   pm2 = m2(w,sigma)
522
   dpm2 = derm2(w,sigma)
523
   VaR = - location - qnorm(alpha)*sqrt(pm2)
524
   derVaR = - as.vector(mu)- qnorm(alpha)*(0.5*as.vector(dpm2))/sqrt(pm2);
525
   contrib = derVaR*as.vector(w) 
526
   return(list( round( VaR , precision ) , round ( contrib , precision ) , round( contrib/VaR , precision) )) 
527
}
528

529
PortgausES =  function(alpha,w,mu,sigma,precision=4){
530
   location = t(w)%*%mu
531
   pm2 = m2(w,sigma)
532
   dpm2 = derm2(w,sigma)
533
   ES = - location + dnorm(qnorm(alpha))*sqrt(pm2)/alpha
534
   derES = - as.vector(mu) + (1/alpha)*dnorm(qnorm(alpha))*(0.5*as.vector(dpm2))/sqrt(pm2);
535
   contrib = as.vector(w)*derES;
536
   return(list( round( ES , precision ) , round( contrib , precision) , round( contrib/ES , precision) )) 
537
}
538

539
PortSkew =  function(w,sigma,M3)
540
{
541
   pm2 = m2(w,sigma)
542
   pm3 = m3(w,M3)
543
   skew = pm3 / pm2^(3/2);
544
   return( skew )
545
}
546

547
PortKurt =  function(w,sigma,M4)
548
{
549
   pm2 = m2(w,sigma)
550
   pm4 = m4(w,M4)
551
   kurt = pm4 / pm2^(2) ;
552
   return( kurt )
553
}
554

555
PortMVaR =  function(alpha,w,mu,sigma,M3,M4,precision=4)
556
{
557
   z = qnorm(alpha)
558
   location = t(w)%*%mu
559
   pm2 = m2(w,sigma)
560
   dpm2 = as.vector( derm2(w,sigma) )
561
   pm3 = m3(w,M3)
562
   dpm3 = as.vector( derm3(w,M3) )
563
   pm4 = m4(w,M4)
564
   dpm4 = as.vector( derm4(w,M4) )
565

566
   skew = pm3 / pm2^(3/2);
567
   exkurt = pm4 / pm2^(2) - 3;
568

569
   derskew = ( 2*(pm2^(3/2))*dpm3 - 3*pm3*sqrt(pm2)*dpm2 )/(2*pm2^3)
570
   derexkurt = ( (pm2)*dpm4 - 2*pm4*dpm2    )/(pm2^3)
571

572
   h = z + (1/6)*(z^2 -1)*skew 
573
   h = h + (1/24)*(z^3 - 3*z)*exkurt - (1/36)*(2*z^3 - 5*z)*skew^2; 
574

575
   MVaR = - location - h*sqrt(pm2)
576

577
   derGausVaR = - as.vector(mu)- qnorm(alpha)*(0.5*as.vector(dpm2))/sqrt(pm2);
578
   derMVaR = derGausVaR + (0.5*dpm2/sqrt(pm2))*( -(1/6)*(z^2 -1)*skew  - (1/24)*(z^3 - 3*z)*exkurt + (1/36)*(2*z^3 - 5*z)*skew^2 )
579
   derMVaR = derMVaR + sqrt(pm2)*( -(1/6)*(z^2 -1)*derskew  - (1/24)*(z^3 - 3*z)*derexkurt + (1/36)*(2*z^3 - 5*z)*2*skew*derskew  )
580
   contrib = as.vector(w)*as.vector(derMVaR)
581
   return(list(  round( MVaR , precision) , round( contrib , precision ), round (contrib/MVaR , precision ) ) ) 
582
}
583

584
derIpower = function(power,h)
585
{
586

587
   fullprod = 1;
588

589
   if( (power%%2)==0 ) #even number: number mod is zero
590
   {
591
      pstar = power/2;
592
      for(j in c(1:pstar))
593
      {
594
         fullprod = fullprod*(2*j)
595
      }
596
      I = -fullprod*h*dnorm(h);
597

598
      for(i in c(1:pstar) )
599
      { 
600
         prod = 1;
601
         for(j in c(1:i) )
602
         {
603
            prod = prod*(2*j)
604
         }
605
         I = I + (fullprod/prod)*(h^(2*i-1))*(2*i-h^2)*dnorm(h)
606
      }
607
   }else{
608
      pstar = (power-1)/2
609
      for(j in c(0:pstar) )
610
      {
611
         fullprod = fullprod*( (2*j)+1 )
612
      }
613
      I = -fullprod*dnorm(h);
614

615
      for(i in c(0:pstar) )
616
      { 
617
         prod = 1;
618
         for(j in c(0:i) )
619
         {
620
            prod = prod*( (2*j) + 1 )
621
         }
622
         I = I + (fullprod/prod)*(h^(2*i)*(2*i+1-h^2) )*dnorm(h)
623
      }
624
   }
625
   return(I)
626
}
627

628

629
PortMES =  function(alpha,w,mu,sigma,M3,M4,precision=4)
630
{
631
   z = qnorm(alpha)
632
   location = t(w)%*%mu
633
   pm2 = m2(w,sigma)
634
   dpm2 = as.vector( derm2(w,sigma) )
635
   pm3 = m3(w,M3)
636
   dpm3 = as.vector( derm3(w,M3) )
637
   pm4 = m4(w,M4)
638
   dpm4 = as.vector( derm4(w,M4) )
639

640
   skew = pm3 / pm2^(3/2);
641
   exkurt = pm4 / pm2^(2) - 3;
642

643
   derskew = ( 2*(pm2^(3/2))*dpm3 - 3*pm3*sqrt(pm2)*dpm2 )/(2*pm2^3)
644
   derexkurt = ( (pm2)*dpm4 - 2*pm4*dpm2    )/(pm2^3)
645

646
   h = z + (1/6)*(z^2 -1)*skew 
647
   h = h + (1/24)*(z^3 - 3*z)*exkurt - (1/36)*(2*z^3 - 5*z)*skew^2; 
648

649
   derh = (1/6)*(z^2 -1)*derskew + (1/24)*(z^3 - 3*z)*derexkurt - (1/18)*(2*z^3 - 5*z)*skew*derskew
650

651
   E = dnorm(h)
652
   E = E + (1/24)*(   Ipower(4,h) - 6*Ipower(2,h) + 3*dnorm(h)   )*exkurt
653
   E = E +  (1/6)*(   Ipower(3,h) - 3*Ipower(1,h)   )*skew;
654
   E = E + (1/72)*(  Ipower(6,h) -15*Ipower(4,h)+ 45*Ipower(2,h) - 15*dnorm(h) )*(skew^2)
655
   E = E/alpha
656
   MES = - location + sqrt(pm2)*E
657

658
   derMES = -mu + 0.5*(dpm2/sqrt(pm2))*E
659
   derE = (1/24)*(   Ipower(4,h) - 6*Ipower(2,h) + 3*dnorm(h)   )*derexkurt 
660
   derE = derE +  (1/6)*(   Ipower(3,h) - 3*Ipower(1,h)   )*derskew
661
   derE = derE + (1/36)*(  Ipower(6,h) -15*Ipower(4,h)+ 45*Ipower(2,h) - 15*dnorm(h) )*skew*derskew
662
   X = -h*dnorm(h) + (1/24)*(  derIpower(4,h) - 6*derIpower(2,h) -3*h*dnorm(h)  )*exkurt 
663
   X = X + (1/6)*( derIpower(3,h) - 3*derIpower(1,h) )*skew 
664
   X = X + (1/72)*( derIpower(6,h) - 15*derIpower(4,h) + 45*derIpower(2,h) + 15*h*dnorm(h)  )*skew^2
665
   derE = derE+derh*X  # X is a scalar
666
   derE = derE/alpha
667
   derMES = derMES + sqrt(pm2)*derE
668
   contrib = as.vector(w)*as.vector(derMES)
669
   return(list(  round( MES , precision ) , round( contrib , precision ), round( contrib/MES, precision )) ) 
670
}
671

672

673
operPortMES =  function(alpha,w,mu,sigma,M3,M4,precision=4)
674
{
675
   z = qnorm(alpha)
676
   location = t(w)%*%mu
677
   pm2 = m2(w,sigma)
678
   dpm2 = as.vector( derm2(w,sigma) )
679
   pm3 = m3(w,M3)
680
   dpm3 = as.vector( derm3(w,M3) )
681
   pm4 = m4(w,M4)
682
   dpm4 = as.vector( derm4(w,M4) )
683

684
   skew = pm3 / pm2^(3/2);
685
   exkurt = pm4 / pm2^(2) - 3;
686

687
   derskew = ( 2*(pm2^(3/2))*dpm3 - 3*pm3*sqrt(pm2)*dpm2 )/(2*pm2^3)
688
   derexkurt = ( (pm2)*dpm4 - 2*pm4*dpm2    )/(pm2^3)
689

690
   h = z + (1/6)*(z^2 -1)*skew 
691
   h = h + (1/24)*(z^3 - 3*z)*exkurt - (1/36)*(2*z^3 - 5*z)*skew^2; 
692
   I1 = Ipower(1,h); I2 = Ipower(2,h); I3 = Ipower(3,h); I4 = Ipower(4,h);  I6 = Ipower(6,h); 
693

694
   derh = (1/6)*(z^2 -1)*derskew + (1/24)*(z^3 - 3*z)*derexkurt - (1/18)*(2*z^3 - 5*z)*skew*derskew
695

696
   E = dnorm(h)
697
   E = E + (1/24)*(   I4 - 6*I2 + 3*dnorm(h)   )*exkurt
698
   E = E +  (1/6)*(   I3 - 3*I1   )*skew;
699
   E = E + (1/72)*( I6 -15*I4+ 45*I2 - 15*dnorm(h) )*(skew^2)
700
   E = E/alpha
701

702
   MES = - location - sqrt(pm2)*min(-E,h)
703

704
   if(-E<=h){
705
      derMES = -mu + 0.5*(dpm2/sqrt(pm2))*E
706
      derE = (1/24)*(   I4 - 6*I2 + 3*dnorm(h)   )*derexkurt 
707
      derE = derE +  (1/6)*(   I3 - 3*I1   )*derskew
708
      derE = derE + (1/36)*(  I6 -15*I4 + 45*I2 - 15*dnorm(h) )*skew*derskew
709
      X = -h*dnorm(h) + (1/24)*(  derIpower(4,h) - 6*derIpower(2,h) -3*h*dnorm(h)  )*exkurt 
710
      X = X + (1/6)*( derIpower(3,h) - 3*derIpower(1,h) )*skew 
711
      X = X + (1/72)*( derIpower(6,h) - 15*derIpower(4,h) + 45*derIpower(2,h) + 15*h*dnorm(h)  )*skew^2
712
      derE = derE+derh*X  # X is a scalar
713
      derE = derE/alpha
714
      derMES = derMES + sqrt(pm2)*derE }else{
715
      derMES = -mu - 0.5*(dpm2/sqrt(pm2))*h - sqrt(pm2)*derh ;  }
716
   contrib = as.vector(w)*as.vector(derMES)
717
   return(list( round( MES, precision) , round( contrib , precision ) , round(contrib/MES,precision) ) ) 
718
}
719

720

721
centeredmoment = function(series,power)
722
{
723
   location = mean(series);
724
   out = sum( (series-location)^power  )/length(series);
725
   return(out);
726
}
727

728
operMES =  function(series,alpha,r)
729
{
730
   z = qnorm(alpha)
731
   location = mean(series);
732
   m2 = centeredmoment(series,2)
733
   m3 = centeredmoment(series,3)
734
   m4 = centeredmoment(series,4)
735
   skew = m3 / m2^(3/2);
736
   exkurt = m4 / m2^(2) - 3;
737

738
   h = z + (1/6)*(z^2 -1)*skew 
739
   if(r==2){ h = h + (1/24)*(z^3 - 3*z)*exkurt - (1/36)*(2*z^3 - 5*z)*skew^2};
740

741
   MES = dnorm(h)
742
   MES = MES + (1/24)*(   Ipower(4,h) - 6*Ipower(2,h) + 3*dnorm(h)   )*exkurt
743
   MES = MES +  (1/6)*(   Ipower(3,h) - 3*Ipower(1,h)   )*skew;
744
   MES = MES + (1/72)*(  Ipower(6,h) -15*Ipower(4,h)+ 45*Ipower(2,h) - 15*dnorm(h) )*(skew^2)
745
   MES = - location - (sqrt(m2))*min( -MES/alpha , h )
746
   return(MES)
747
}
748

749
TwoVarPlot <- function(xvar, y1var, y2var, labels, noincs = 5,marks=c(1,2), legpos, leglabs, title)
750
{ 
751
   # https://stat.ethz.ch/pipermail/r-help/2000-September/008182.html
752

753
   # plots an x y1 y2 using left and right axes for the different y's
754
   # rescales y2 to fit in the same space as y1
755

756
   # noincs - no of divisions in the axis labels
757
   # marks - type of marker for each y
758
   # legpos - legend position
759
   # leglabs - legend labels
760
  
761
   # rescale to fit on same axis
762
   scaledy2var <- (y2var - min(y2var)) / (max(y2var) - min(y2var))
763
   scaledy2var <- (scaledy2var * (max(y1var) - min(y1var))) + min(y1var)
764

765
   # plot it up and add the points
766
   plot(xvar, y1var, xlab=labels[1], ylab="", axes=F, pch=marks[1],main=title,type="l")
767
   lines(xvar, scaledy2var, lty=3 )
768

769
   # make up some labels and positions
770
   y1labs <- round(seq(min(y1var), max(y1var), length=noincs),2)
771

772
   # convert these to the y2 axis scaling
773
   y2labs <- (y1labs - min(y1var)) / (max(y1var) - min(y1var))
774
   y2labs <- (y2labs * (max(y2var) - min(y2var))) + min(y2var)
775
   y2labs <- round(y2labs, 2)
776

777
   axis(1)
778
   axis(2, at=y1labs, labels=y1labs)
779
   axis(4, at=y1labs, labels=y2labs)
780
   mtext(labels[3], side=4, line=2)
781
   mtext(labels[2], side=2, line=2)
782
   box()
783

784
   legend( legend=leglabs, lty = c(1,3), bty="o", x=legpos)
785
}
786

787

788

789
Product

Resources

Company
