CoCalc -- Risk_budget

GitHub Repository: braverock/portfolioanalytics
Path: blob/master/sandbox/paper_analysis/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);  
19
     require( DEoptim )
20

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

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

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

51
    if( is.null(penalty) ){ 
52
       abspenalty = 1e6; relpenalty = 100*N 
53
    }else{
54
       abspenalty = relpenalty = penalty; 
55
    }; 
56

57
    if( optimize_method == "quadprog" & minriskcriterion == "StdDev"  ){
58
       require(quadprog)
59
       dvec = matrix( rep(0,N) , ncol=1 )
60
       # only risky assets
61
       sumweights = 1;
62
       Amat = rep(1,N);
63
       bvec = sumweights;
64
       # No short selling
65
       minweight = 0 ;
66
       # first round (afterwards, we impose minassetweight on assets with non-zero weight [heuristic]
67
       Amat = rbind( Amat ,  diag(x =1,nrow=N,ncol=N) ,  diag(x =-1,nrow=N,ncol=N));
68
       # you need to allow for an epsilon to make it feasible in the second stage
69
       bvec = matrix( c( bvec , rep(minweight,N), -upper  ),ncol=1)
70
       minw = solve.QP( Dmat = sigma   , dvec = dvec , Amat = t(Amat) , bvec = bvec  , meq =1  )$solution
71
    } 
72

73
    if( optimize_method == "DEoptim" | optimize_method == "DEoptim+L-BFGS-B" ){
74
        objective = function( w ){
75
           if( sum(w) < 1e-4 ){
76
               w <- w + 1e-2
77
           }
78
           w = w/sum(w) 
79
           cont = percriskcontrib( w );  percrisk = cont[[3]]; crisk = cont[[2]] ; 
80
           if(MinMaxComp){ out = max( crisk ) }else{ out = prisk(w) }
81
           # add weight constraints
82
           out_con = out;
83
           if( any(lower !=0) | any(upper!=1) ){
84
              out_con = out_con + out*abspenalty*sum( 1*( w < lower )+1*( w > upper ) ) # can happen because of the standardization
85
           }
86
           # add risk budget constraint 
87
           if( any(RBupper < Inf ) ){
88
               con1 = sum(  (percrisk-RBupper)[RBupper < Inf]*( percrisk > RBupper  )[RBupper < Inf] ,na.rm=TRUE ) ; 
89
               out_con = out_con + out*relpenalty*con1
90
           }
91
           if( any(RBlower > -Inf ) ){
92
               con2 = sum(  (RBlower-percrisk)[RBlower > -Inf]*( percrisk < RBlower  )[RBlower > -Inf] ,na.rm=TRUE ) ; 
93
               out_con = out_con + out*relpenalty*con2
94
           }
95
           # add minimum risk constraint
96
           if( Riskupper < Inf ){
97
              con = ( prisk(w) - Riskupper)*( prisk(w) > Riskupper ); 
98
              out_con = out_con + out*relpenalty*con
99
           }
100
           # add minimum return constraint 
101
           # portfolio return and risk are in the same unit, but portfolio return is typically some orders of magnitude smaller
102
           # say: as a very conservative choice: 100 
103
           preturn = sum( w*mu ) ;
104
           if( Returnlower > -Inf ){
105
              con = ( Returnlower - preturn)*( preturn < Returnlower ); if( is.na(con) ){ con = 0 }
106
              out_con = out_con + out*100*relpenalty*con
107
           }
108
           #der = derprisk(w)[w>0.01]
109
           #der = abs(der -mean( der ) )
110
           #if( sum( der )>0.01 ){
111
           #   out_con = out_con + out/sum( der )* max( der )}
112
           return(out_con)
113
         }
114
         set.seed(1234)
115
         if( is.null(controlDE) ){
116
            # initial population is generated with random_portfolios function
117
            require('PortfolioAnalytics')
118
            eps = 0.01
119
            rpconstraint<-constraint(assets=length(lower), min_sum=(1-eps), max_sum=(1+eps), 
120
              min=lower, max=upper, weight_seq=generatesequence())
121
            rp<- random_portfolios(rpconstraints=rpconstraint,permutations=200)
122
            controlDE <- list(reltol=1e-6,steptol=150, itermax = 5000,trace = 100,
123
                strategy=6, c=.4,
124
                NP=as.numeric(nrow(rp)),initialpop=rp) 
125
         }
126
         # print(controlDE)
127
         minw = DEoptim( fn = objective ,  lower = lower , upper = upper , control = controlDE)
128
         fvalues = minw$member$bestval
129
         #diff = as.vector( quantile(fvalues,0.10) - min(fvalues) )
130
         #print(c("diff",diff))
131
         best = min( fvalues ) ; print(best)
132
         bestsol = minw ; 
133
         minw = bestsol$optim$bestmem/sum(bestsol$optim$bestmem)  ; #full investment constraint
134
         print( bestsol$optim$iter )
135
    }
136
   
137

138
    if( optimize_method == "DEoptim+L-BFGS-B" ){
139
         print(minw) ; print( objective(minw) ) ;         
140
         print("local improvement")
141
         out = optim( par=minw, f = objective, lower = lower, upper = upper, method="L-BFGS-B" ) 
142
         minw = out$par/sum(out$par)
143
         print(minw) ; print( objective(minw) )
144
    }
145

146
    if( optimize_method == "L-BFGS-B" ){
147
        objective = function( w ){
148
           if(sum(w)==0){w=w+0.001}
149
           w = w/sum(w)  
150
           cont = percriskcontrib( w );  percrisk = cont[[3]]; crisk = cont[[2]] ; 
151
           if(MinMaxComp){ out = max( crisk ) }else{ out = prisk(w) }
152
           # add risk budget constraint 
153
           if( max(RBupper) != Inf ){ con1 = sum(  (percrisk-RBupper)*( percrisk > RBupper  ),na.rm=TRUE )}else{ con1 = 0 } 
154
           if( min(RBlower) != -Inf ){con2 = sum(  (RBlower-percrisk)*( percrisk < RBlower  ),na.rm=TRUE )}else{ con2 = 0 } 
155
           con = con1 + con2; if( is.na(con) ){ con = 0 }  
156
           out_con = out + out*relpenalty*con
157
           # add minimum risk constraint
158
           if( Riskupper != Inf ){ con = ( prisk(w) - Riskupper)*( prisk(w) > Riskupper )}else{ con = 0 }; 
159
           if( is.na(con) ){ con = 0 }  
160
           out_con = out_con + out*relpenalty*con
161
           return(out_con)
162
         }
163

164
         if(Returnlower==-Inf){ 
165
            inittheta = rep(1,N)/N 
166
            out = optim( par=inittheta, f = objective, lower = lower, upper = upper ) 
167
         }else{
168
            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))
169
            inittheta = rep(0.001,N);
170
            inittheta[mu==max(mu)] = 1; inittheta = 1-sum(inittheta[mu!=max(mu)]  );   
171
            out = constrOptim( theta=inittheta, f = objective, grad=NULL,ui=Amat,
172
                            ci = c(rep(0,N),rep(-1,N),0.99999,-1.0001,Returnlower) ) 
173
         }
174

175
         minw = out$par/sum(out$par)
176
    }
177
    cont = percriskcontrib( minw );  percrisk = cont[[3]];    crisk = cont[[2]] ;    
178
    # check 
179
    print( "out = list( minw , sum( minw*mu ) , prisk(minw) , percriskcontrib(minw)" )
180
    out = list( minw , sum( minw*mu ) , prisk(minw) , percrisk , crisk )
181
    print( out )
182
    return(out)
183
}
184

185
findportfolio.dynamic = function(R, mulist , sigmalist, M3list , M4list , from, to, names.input = NA,  names.assets, 
186
              	p = 0.95 , priskbudget = 0.95,  mincriterion = "mES" , percriskcontribcriterion = "mES"  ,
187
                  strategy , controlDE =  NULL , optimize_method = "DEoptim+L-BFGS-B" )
188
{ # @author Kris Boudt and Brian G. Peterson
189

190
    # Description:
191
    #
192
    # Performs a loop over the reallocation periods with estimation samples given by from:to
193
    # It calls the function RBconportfolio to obtain the optimal weights of the strategy.
194
    #
195
    # @todo
196
    #
197
    # R                 matrix/zoo holding historical returns on risky assets
198
    #
199
    # names             vector holding the names of the .csv files to be read
200
    #
201
    # from, to          define the estimation sample
202
    #
203
    # criteria	      the criterion to be optimized
204
    #
205
    # columns.crit      the columns of R in which the criteria are located
206
    #
207
    # percriskcontribcriterion risk measure used for the risk budget constraints
208
    #
209
    # strategy = c( "EqualRisk" , "EqualWeight" , "MinRisk" , "MinRiskConc" , 
210
    #                "MinRisk_PositionLimit"     , "MinRisk_RiskLimit"     , "MinRisk_ReturnTarget",
211
    #                "MinRiskConc_PositionLimit" , "MinRiskConc_RiskLimit" , "MinRiskConc_ReturnTarget")
212

213
    # Return:
214
    # List with first element optimal weights per reallocation period and associated percentage CVaR contributions. 
215

216
    # Create a matrix that will hold for each method and each vector the best weights
217

218
    cPeriods = length(from);    
219

220
    out = matrix(  rep(0, cPeriods*(cAssets)) , ncol= (cAssets) );
221
    RCout = matrix(  rep(0, cPeriods*(cAssets)) , ncol= (cAssets) );
222
    RiskReturnout = matrix(  rep(0, cPeriods*2 ) , ncol= 2 );
223
 
224
    # first cPeriods rows correspond to cCriteria[1] and so on
225

226
    # downside risk
227
    alpha = 1 - p;
228
    alphariskbudget = 1-priskbudget;
229

230
    # Estimation of the return mean vector, covariance, coskewness and cokurtosis matrix
231

232

233
    if(strategy=="EqualRisk"){
234
       lower = rep(0,cAssets); upper=rep(1,cAssets) 
235
       RBlower = rep(1/cAssets,cAssets) ; RBupper = rep(1/cAssets,cAssets)  ;
236
    }
237

238
    if(strategy=="EqualWeight"){
239
       lower = rep(1/cAssets,cAssets); upper=rep(1/cAssets,cAssets) 
240
       RBlower = rep(-Inf,cAssets) ; RBupper = rep(Inf,cAssets)  ;
241
    }
242

243
    if(strategy=="MinRisk" | strategy=="MinRiskConc" | strategy=="MinRisk_ReturnTarget" | strategy=="MinRiskConc_ReturnTarget"){
244
       lower = rep(0,cAssets); upper=rep(1,cAssets) 
245
       RBlower = rep(-Inf,cAssets) ; RBupper = rep(Inf,cAssets)  ;
246
    }
247

248
    MinMaxComp = F; mutarget = -Inf; 
249
    if( strategy=="MinRiskConc" | strategy=="MinRiskConc_PositionLimit" | strategy=="MinRiskConc_RiskLimit" | strategy=="MinRiskConc_ReturnTarget"  ){
250
        MinMaxComp = T;
251
    }
252

253
   if(strategy=="MinRisk_PositionLimit" | strategy=="MinRiskConc_PositionLimit"){
254
      lower = rep(0,cAssets); upper=rep(0.4,cAssets) 
255
      RBlower = rep(-Inf,cAssets) ; RBupper = rep(Inf,cAssets)  ;
256
   }
257

258
   if(strategy=="MinRisk_RiskLimit" | strategy=="MinRiskConc_RiskLimit"){
259
      lower = rep(0,cAssets); upper=rep(1,cAssets) 
260
      RBlower = rep(-Inf,cAssets) ; RBupper = rep(0.40,cAssets)  ;
261
   }
262

263
   # initial population is generated with random_portfolios function
264
   if(is.null(controlDE) & optimize_method != "quadprog" ){
265
     require('PortfolioAnalytics')
266
     eps = 0.025
267
     rpconstraint<-constraint(assets=length(lower), min_sum=(1-eps), max_sum=(1+eps), 
268
             min=lower, max=upper, weight_seq=generatesequence())
269
     rp<- random_portfolios(rpconstraints=rpconstraint,permutations=controlDE$NP)
270
     controlDE = list(reltol=1e-6,steptol=150, itermax = 5000,trace = 500, strategy=6, c=.4,
271
                NP=as.numeric(nrow(rp)),initialpop=rp) 
272
   }
273
    for( per in c(1:cPeriods) ){
274

275
       print("-----------New estimation period ends on observation------------------")
276
       print( paste(to[per],"out of total number of obs equal to", max(to) ));
277
       print("----------------------------------------------------------------")
278

279
       mu    = 100*mulist[[per]] ; 
280
       sigma = (100^2)*sigmalist[[per]] ; 
281
       M3    = (100^3)*M3list[[per]] ; 
282
       M4    = (100^4)*M4list[[per]] ; 
283

284
     
285
       mESfun = function(series){ return( operMES(series,alpha=alpha,2) ) }
286
       
287
 
288
       if(strategy=="MinRisk_ReturnTarget" | strategy=="MinRiskConc_ReturnTarget"){
289
           mutarget = mean( mu );
290
           print( c("Minimum return requirement is" , mutarget) )
291
       }
292

293
       if(strategy=="EqualWeight"){
294
           sol1 = rep(1/cAssets,cAssets);
295
           switch( percriskcontribcriterion ,
296
                StdDev = { percriskcontrib = function(w){ cont = Portsd(w,sigma=sigma) ; return( cont  ) }},
297
                GVaR = {percriskcontrib = function(w){ cont =  PortgausVaR(w,alpha=alphariskbudget,mu=mu,sigma=sigma) ; return( cont ) }},
298
                GES = {percriskcontrib = function(w){ cont =  PortgausES(w,mu=mu,alpha=alphariskbudget,sigma=sigma) ; return( cont  ) }},
299
                mVaR = {percriskcontrib = function(w){ cont = PortMVaR(w,mu=mu,alpha=alphariskbudget,sigma=sigma,M3=M3,M4=M4) ; return( cont  ) }},
300
                mES = {percriskcontrib = function(w){ cont = operPortMES(w,mu=mu,alpha=alphariskbudget,sigma=sigma,M3=M3,M4=M4) ; return( cont ) }
301
                }
302
           )
303
           sol2 = percriskcontrib( sol1 )
304
           solution = list( sol1 , sol2[[3]] , sol2[[1]] , mean(mu) ) ;   
305
       }else{
306
           if( per > 1 & optimize_method != "quadprog" ){
307
              controlDE$initialpop = rbind( rp , out[1:(per-1),] )
308
              controlDE$NP         = as.numeric( nrow(controlDE$initialpop) )
309
           }
310
           #controlDE = list(reltol=1e-6,steptol=150, itermax = 5000,trace = 500, strategy=6, c=.4,
311
           #     NP=as.numeric(nrow(rp)),initialpop=rp) 
312
           solution = PortfolioOptim(  minriskcriterion = mincriterion, MinMaxComp = MinMaxComp, 
313
                                       percriskcontribcriterion = percriskcontribcriterion , 
314
                                       mu = mu , sigma = sigma, M3=M3 , M4=M4 , alpha = alpha , alphariskbudget = alphariskbudget , 
315
                                       lower = lower , upper = upper , Riskupper = Inf , Returnlower= mutarget , RBlower = RBlower, RBupper = RBupper , 
316
                                       controlDE = controlDE , optimize_method = optimize_method ) 
317
           # [[1]] : weights ; [[2]] : expected portfolio return ; [[3]] : portfolio CVaR
318
           # [[4]] : percentage CVaR ; [[5]] : component CVaR
319
           solution = list( solution[[1]] , solution[[4]] , solution[[3]]/100 , solution[[2]]/100 ); 
320
       }
321
       out[ per, ]    = as.vector( solution[[1]] ) 
322
       RCout[per,  ]  = as.vector( solution[[2]] )   
323
       RiskReturnout[per,] = c( as.numeric(solution[[3]]) , as.numeric(solution[[4]]) )
324

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

327
  # Output save
328
  rownames(out) = rownames(RCout) = rownames(RiskReturnout) =  names.input; 
329
  colnames(out) = colnames(RCout) = names.assets;
330
  colnames(RiskReturnout) = c( "risk" , "expreturn" )
331

332
  EWweights = c( rep(1/cAssets,cAssets) )
333
  EWweights = matrix ( rep(EWweights,cPeriods) , ncol=(cAssets) , byrow = TRUE )
334
  rownames(EWweights) = names.input; colnames(EWweights) = names.assets;
335

336
  return( list( out, RCout, RiskReturnout) ) 
337
}
338

339
clean.boudt2 = 
340
function (R, alpha = 0.01, trim = 0.001) 
341
{
342
    # @author Kris Boudt and Brian Peterson
343
    # Cleaning method as described in 
344
    #     Boudt, Peterson and Croux. 2009. Estimation and decomposition of downside risk for portfolios with non-normal returns. Journal of Risk.
345

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

380

381
Ipower = function(power,h)
382
{
383

384
   fullprod = 1;
385

386
   if( (power%%2)==0 ) #even number: number mod is zero
387
   {
388
      pstar = power/2;
389
      for(j in c(1:pstar))
390
      {
391
         fullprod = fullprod*(2*j)
392
      }
393
      I = fullprod*dnorm(h);
394

395
      for(i in c(1:pstar) )
396
      { 
397
         prod = 1;
398
         for(j in c(1:i) )
399
         {
400
            prod = prod*(2*j)
401
         }
402
         I = I + (fullprod/prod)*(h^(2*i))*dnorm(h)
403
      }
404
   }
405
   else{
406
      pstar = (power-1)/2
407
      for(j in c(0:pstar) )
408
      {
409
         fullprod = fullprod*( (2*j)+1 )
410
      }
411
      I = -fullprod*pnorm(h);
412

413
      for(i in c(0:pstar) )
414
      { 
415
         prod = 1;
416
         for(j in c(0:i) )
417
         {
418
            prod = prod*( (2*j) + 1 )
419
         }
420
         I = I + (fullprod/prod)*(h^(  (2*i) + 1))*dnorm(h)
421
      }
422
   }
423
   return(I)
424
}
425

426

427
# Definition of statistics needed to compute Gaussian and modified VaR and ES for the return series of portfolios
428
# and to compute the contributions to portfolio downside risk, made by the different positions in the portfolio.
429
#----------------------------------------------------------------------------------------------------------------
430

431

432

433
m2 = function(w,sigma)
434
{
435
   return(t(w)%*%sigma%*%w)
436
}
437
derm2 = function(w,sigma)
438
{
439
   return(2*sigma%*%w)
440
}
441
m3 = function(w,M3)
442
{
443
   return(t(w)%*%M3%*%(w%x%w))
444
}
445
derm3 = function(w,M3)
446
{
447
   return(3*M3%*%(w%x%w))
448
}
449
m4 = function(w,M4)
450
{
451
   return(t(w)%*%M4%*%(w%x%w%x%w))
452
}
453
derm4 = function(w,M4)
454
{
455
   return(4*M4%*%(w%x%w%x%w))
456
}
457

458
StdDevfun = function(w,sigma){ return(  sqrt( t(w)%*%sigma%*%w  )) }
459

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

462
mVaRfun = function(w,alpha,mu,sigma,M3,M4){
463
    pm4 = t(w)%*%M4%*%(w%x%w%x%w) ; pm3 = t(w)%*%M3%*%(w%x%w) ; pm2 =  t(w)%*%sigma%*%w ; 
464
    skew = pm3 / pm2^(3/2);
465
    exkurt = pm4 / pm2^(2) - 3; z = qnorm(alpha);
466
    h = z + (1/6)*(z^2 -1)*skew
467
    h = h + (1/24)*(z^3 - 3*z)*exkurt - (1/36)*(2*z^3 - 5*z)*skew^2;
468
    return (- (t(w)%*%mu) - h*sqrt( pm2  ) ) }
469

470
resmVaRfun = function(w,alpha,mu,sigma,ressigma,M3,M4){
471
    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 ;
472
    skew = pm3 / respm2^(3/2);
473
    exkurt = pm4 / respm2^(2) - 3; z = qnorm(alpha);
474
    h = z + (1/6)*(z^2 -1)*skew
475
    h = h + (1/24)*(z^3 - 3*z)*exkurt - (1/36)*(2*z^3 - 5*z)*skew^2;
476
    return (- (t(w)%*%mu) - h*sqrt( pm2  ) ) }
477

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

481
operMESfun = function(w,alpha,mu,sigma,M3,M4){
482
    pm4 = t(w)%*%M4%*%(w%x%w%x%w) ; pm3 = t(w)%*%M3%*%(w%x%w) ; pm2 =  t(w)%*%sigma%*%w ; 
483
    skew = pm3 / pm2^(3/2);
484
    exkurt = pm4 / pm2^(2) - 3; z = qnorm(alpha);
485
    h = z + (1/6)*(z^2 -1)*skew
486
    h = h + (1/24)*(z^3 - 3*z)*exkurt - (1/36)*(2*z^3 - 5*z)*skew^2;
487
    E = dnorm(h)
488
    E = E + (1/24)*(   Ipower(4,h) - 6*Ipower(2,h) + 3*dnorm(h)   )*exkurt
489
    E = E +  (1/6)*(   Ipower(3,h) - 3*Ipower(1,h)   )*skew;
490
    E = E + (1/72)*(  Ipower(6,h) -15*Ipower(4,h)+ 45*Ipower(2,h) - 15*dnorm(h) )*(skew^2)
491
    E = E/alpha
492
    return (- (t(w)%*%mu) - sqrt(pm2)*min(-E,h) )
493
}
494

495

496
Portmean = function(w,mu,precision=4)
497
{
498
   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) )
499
}
500

501
Portsd =  function(w,sigma,precision=4)
502
{
503
   pm2 = m2(w,sigma)
504
   dpm2 = derm2(w,sigma)
505
   dersd = (0.5*as.vector(dpm2))/sqrt(pm2);
506
   contrib = dersd*as.vector(w)
507
   return(list(  round( sqrt(pm2) , precision ) , round( contrib , precision ) , round ( contrib/sqrt(pm2) , precision) )) 
508
}
509

510

511
PortgausVaR =  function(alpha,w,mu,sigma,precision=4){
512
   location = t(w)%*%mu
513
   pm2 = m2(w,sigma)
514
   dpm2 = derm2(w,sigma)
515
   VaR = - location - qnorm(alpha)*sqrt(pm2)
516
   derVaR = - as.vector(mu)- qnorm(alpha)*(0.5*as.vector(dpm2))/sqrt(pm2);
517
   contrib = derVaR*as.vector(w) 
518
   return(list( round( VaR , precision ) , round ( contrib , precision ) , round( contrib/VaR , precision) )) 
519
}
520

521
PortgausES =  function(alpha,w,mu,sigma,precision=4){
522
   location = t(w)%*%mu
523
   pm2 = m2(w,sigma)
524
   dpm2 = derm2(w,sigma)
525
   ES = - location + dnorm(qnorm(alpha))*sqrt(pm2)/alpha
526
   derES = - as.vector(mu) + (1/alpha)*dnorm(qnorm(alpha))*(0.5*as.vector(dpm2))/sqrt(pm2);
527
   contrib = as.vector(w)*derES;
528
   return(list( round( ES , precision ) , round( contrib , precision) , round( contrib/ES , precision) )) 
529
}
530

531
PortSkew =  function(w,sigma,M3)
532
{
533
   pm2 = m2(w,sigma)
534
   pm3 = m3(w,M3)
535
   skew = pm3 / pm2^(3/2);
536
   return( skew )
537
}
538

539
PortKurt =  function(w,sigma,M4)
540
{
541
   pm2 = m2(w,sigma)
542
   pm4 = m4(w,M4)
543
   kurt = pm4 / pm2^(2) ;
544
   return( kurt )
545
}
546

547
PortMVaR =  function(alpha,w,mu,sigma,M3,M4,precision=4)
548
{
549
   z = qnorm(alpha)
550
   location = t(w)%*%mu
551
   pm2 = m2(w,sigma)
552
   dpm2 = as.vector( derm2(w,sigma) )
553
   pm3 = m3(w,M3)
554
   dpm3 = as.vector( derm3(w,M3) )
555
   pm4 = m4(w,M4)
556
   dpm4 = as.vector( derm4(w,M4) )
557

558
   skew = pm3 / pm2^(3/2);
559
   exkurt = pm4 / pm2^(2) - 3;
560

561
   derskew = ( 2*(pm2^(3/2))*dpm3 - 3*pm3*sqrt(pm2)*dpm2 )/(2*pm2^3)
562
   derexkurt = ( (pm2)*dpm4 - 2*pm4*dpm2    )/(pm2^3)
563

564
   h = z + (1/6)*(z^2 -1)*skew 
565
   h = h + (1/24)*(z^3 - 3*z)*exkurt - (1/36)*(2*z^3 - 5*z)*skew^2; 
566

567
   MVaR = - location - h*sqrt(pm2)
568

569
   derGausVaR = - as.vector(mu)- qnorm(alpha)*(0.5*as.vector(dpm2))/sqrt(pm2);
570
   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 )
571
   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  )
572
   contrib = as.vector(w)*as.vector(derMVaR)
573
   return(list(  round( MVaR , precision) , round( contrib , precision ), round (contrib/MVaR , precision ) ) ) 
574
}
575

576
derIpower = function(power,h)
577
{
578

579
   fullprod = 1;
580

581
   if( (power%%2)==0 ) #even number: number mod is zero
582
   {
583
      pstar = power/2;
584
      for(j in c(1:pstar))
585
      {
586
         fullprod = fullprod*(2*j)
587
      }
588
      I = -fullprod*h*dnorm(h);
589

590
      for(i in c(1:pstar) )
591
      { 
592
         prod = 1;
593
         for(j in c(1:i) )
594
         {
595
            prod = prod*(2*j)
596
         }
597
         I = I + (fullprod/prod)*(h^(2*i-1))*(2*i-h^2)*dnorm(h)
598
      }
599
   }else{
600
      pstar = (power-1)/2
601
      for(j in c(0:pstar) )
602
      {
603
         fullprod = fullprod*( (2*j)+1 )
604
      }
605
      I = -fullprod*dnorm(h);
606

607
      for(i in c(0:pstar) )
608
      { 
609
         prod = 1;
610
         for(j in c(0:i) )
611
         {
612
            prod = prod*( (2*j) + 1 )
613
         }
614
         I = I + (fullprod/prod)*(h^(2*i)*(2*i+1-h^2) )*dnorm(h)
615
      }
616
   }
617
   return(I)
618
}
619

620

621
PortMES =  function(alpha,w,mu,sigma,M3,M4,precision=4)
622
{
623
   z = qnorm(alpha)
624
   location = t(w)%*%mu
625
   pm2 = m2(w,sigma)
626
   dpm2 = as.vector( derm2(w,sigma) )
627
   pm3 = m3(w,M3)
628
   dpm3 = as.vector( derm3(w,M3) )
629
   pm4 = m4(w,M4)
630
   dpm4 = as.vector( derm4(w,M4) )
631

632
   skew = pm3 / pm2^(3/2);
633
   exkurt = pm4 / pm2^(2) - 3;
634

635
   derskew = ( 2*(pm2^(3/2))*dpm3 - 3*pm3*sqrt(pm2)*dpm2 )/(2*pm2^3)
636
   derexkurt = ( (pm2)*dpm4 - 2*pm4*dpm2    )/(pm2^3)
637

638
   h = z + (1/6)*(z^2 -1)*skew 
639
   h = h + (1/24)*(z^3 - 3*z)*exkurt - (1/36)*(2*z^3 - 5*z)*skew^2; 
640

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

643
   E = dnorm(h)
644
   E = E + (1/24)*(   Ipower(4,h) - 6*Ipower(2,h) + 3*dnorm(h)   )*exkurt
645
   E = E +  (1/6)*(   Ipower(3,h) - 3*Ipower(1,h)   )*skew;
646
   E = E + (1/72)*(  Ipower(6,h) -15*Ipower(4,h)+ 45*Ipower(2,h) - 15*dnorm(h) )*(skew^2)
647
   E = E/alpha
648
   MES = - location + sqrt(pm2)*E
649

650
   derMES = -mu + 0.5*(dpm2/sqrt(pm2))*E
651
   derE = (1/24)*(   Ipower(4,h) - 6*Ipower(2,h) + 3*dnorm(h)   )*derexkurt 
652
   derE = derE +  (1/6)*(   Ipower(3,h) - 3*Ipower(1,h)   )*derskew
653
   derE = derE + (1/36)*(  Ipower(6,h) -15*Ipower(4,h)+ 45*Ipower(2,h) - 15*dnorm(h) )*skew*derskew
654
   X = -h*dnorm(h) + (1/24)*(  derIpower(4,h) - 6*derIpower(2,h) -3*h*dnorm(h)  )*exkurt 
655
   X = X + (1/6)*( derIpower(3,h) - 3*derIpower(1,h) )*skew 
656
   X = X + (1/72)*( derIpower(6,h) - 15*derIpower(4,h) + 45*derIpower(2,h) + 15*h*dnorm(h)  )*skew^2
657
   derE = derE+derh*X  # X is a scalar
658
   derE = derE/alpha
659
   derMES = derMES + sqrt(pm2)*derE
660
   contrib = as.vector(w)*as.vector(derMES)
661
   return(list(  round( MES , precision ) , round( contrib , precision ), round( contrib/MES, precision )) ) 
662
}
663

664
derMES =  function(alpha,w,mu,sigma,M3,M4,precision=4)
665
{
666
   z = qnorm(alpha)
667
   location = t(w)%*%mu
668
   pm2 = m2(w,sigma)
669
   dpm2 = as.vector( derm2(w,sigma) )
670
   pm3 = m3(w,M3)
671
   dpm3 = as.vector( derm3(w,M3) )
672
   pm4 = m4(w,M4)
673
   dpm4 = as.vector( derm4(w,M4) )
674

675
   skew = pm3 / pm2^(3/2);
676
   exkurt = pm4 / pm2^(2) - 3;
677

678
   derskew = ( 2*(pm2^(3/2))*dpm3 - 3*pm3*sqrt(pm2)*dpm2 )/(2*pm2^3)
679
   derexkurt = ( (pm2)*dpm4 - 2*pm4*dpm2    )/(pm2^3)
680

681
   h = z + (1/6)*(z^2 -1)*skew 
682
   h = h + (1/24)*(z^3 - 3*z)*exkurt - (1/36)*(2*z^3 - 5*z)*skew^2; 
683

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

686
   E = dnorm(h)
687
   E = E + (1/24)*(   Ipower(4,h) - 6*Ipower(2,h) + 3*dnorm(h)   )*exkurt
688
   E = E +  (1/6)*(   Ipower(3,h) - 3*Ipower(1,h)   )*skew;
689
   E = E + (1/72)*(  Ipower(6,h) -15*Ipower(4,h)+ 45*Ipower(2,h) - 15*dnorm(h) )*(skew^2)
690
   E = E/alpha
691
   MES = - location + sqrt(pm2)*E
692

693
   derMES = -mu + 0.5*(dpm2/sqrt(pm2))*E
694
   derE = (1/24)*(   Ipower(4,h) - 6*Ipower(2,h) + 3*dnorm(h)   )*derexkurt 
695
   derE = derE +  (1/6)*(   Ipower(3,h) - 3*Ipower(1,h)   )*derskew
696
   derE = derE + (1/36)*(  Ipower(6,h) -15*Ipower(4,h)+ 45*Ipower(2,h) - 15*dnorm(h) )*skew*derskew
697
   X = -h*dnorm(h) + (1/24)*(  derIpower(4,h) - 6*derIpower(2,h) -3*h*dnorm(h)  )*exkurt 
698
   X = X + (1/6)*( derIpower(3,h) - 3*derIpower(1,h) )*skew 
699
   X = X + (1/72)*( derIpower(6,h) - 15*derIpower(4,h) + 45*derIpower(2,h) + 15*h*dnorm(h)  )*skew^2
700
   derE = derE+derh*X  # X is a scalar
701
   derE = derE/alpha
702
   derMES = derMES + sqrt(pm2)*derE
703
   return(derMES)
704
}
705

706

707
operPortMES =  function(alpha,w,mu,sigma,M3,M4)
708
{
709
   z = qnorm(alpha)
710
   location = t(w)%*%mu
711
   pm2 = m2(w,sigma)
712
   dpm2 = as.vector( derm2(w,sigma) )
713
   pm3 = m3(w,M3)
714
   dpm3 = as.vector( derm3(w,M3) )
715
   pm4 = m4(w,M4)
716
   dpm4 = as.vector( derm4(w,M4) )
717

718
   skew = pm3 / pm2^(3/2);
719
   exkurt = pm4 / pm2^(2) - 3;
720

721
   derskew = ( 2*(pm2^(3/2))*dpm3 - 3*pm3*sqrt(pm2)*dpm2 )/(2*pm2^3)
722
   derexkurt = ( (pm2)*dpm4 - 2*pm4*dpm2    )/(pm2^3)
723

724
   h = z + (1/6)*(z^2 -1)*skew 
725
   h = h + (1/24)*(z^3 - 3*z)*exkurt - (1/36)*(2*z^3 - 5*z)*skew^2; 
726
   I1 = Ipower(1,h); I2 = Ipower(2,h); I3 = Ipower(3,h); I4 = Ipower(4,h);  I6 = Ipower(6,h); 
727

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

730
   E = dnorm(h)
731
   E = E + (1/24)*(   I4 - 6*I2 + 3*dnorm(h)   )*exkurt
732
   E = E +  (1/6)*(   I3 - 3*I1   )*skew;
733
   E = E + (1/72)*( I6 -15*I4+ 45*I2 - 15*dnorm(h) )*(skew^2)
734
   E = E/alpha
735

736
   MES = - location - sqrt(pm2)*min(-E,h)
737

738
   if(-E<=h){
739
      derMES = -mu + 0.5*(dpm2/sqrt(pm2))*E
740
      derE = (1/24)*(   I4 - 6*I2 + 3*dnorm(h)   )*derexkurt 
741
      derE = derE +  (1/6)*(   I3 - 3*I1   )*derskew
742
      derE = derE + (1/36)*(  I6 -15*I4 + 45*I2 - 15*dnorm(h) )*skew*derskew
743
      X = -h*dnorm(h) + (1/24)*(  derIpower(4,h) - 6*derIpower(2,h) -3*h*dnorm(h)  )*exkurt 
744
      X = X + (1/6)*( derIpower(3,h) - 3*derIpower(1,h) )*skew 
745
      X = X + (1/72)*( derIpower(6,h) - 15*derIpower(4,h) + 45*derIpower(2,h) + 15*h*dnorm(h)  )*skew^2
746
      derE = derE+derh*X  # X is a scalar
747
      derE = derE/alpha
748
      derMES = derMES + sqrt(pm2)*derE
749
    }else{
750
       derMES = -mu - 0.5*(dpm2/sqrt(pm2))*h - sqrt(pm2)*derh ; 
751
    }
752
    contrib = as.vector(w)*as.vector(derMES)
753
    return(list(  MES ,  contrib , contrib/MES ) ) 
754
}
755

756
deroperPortMES =  function(alpha,w,mu,sigma,M3,M4)
757
{
758
   z = qnorm(alpha)
759
   location = t(w)%*%mu
760
   pm2 = m2(w,sigma)
761
   dpm2 = as.vector( derm2(w,sigma) )
762
   pm3 = m3(w,M3)
763
   dpm3 = as.vector( derm3(w,M3) )
764
   pm4 = m4(w,M4)
765
   dpm4 = as.vector( derm4(w,M4) )
766

767
   skew = pm3 / pm2^(3/2);
768
   exkurt = pm4 / pm2^(2) - 3;
769

770
   derskew = ( 2*(pm2^(3/2))*dpm3 - 3*pm3*sqrt(pm2)*dpm2 )/(2*pm2^3)
771
   derexkurt = ( (pm2)*dpm4 - 2*pm4*dpm2    )/(pm2^3)
772

773
   h = z + (1/6)*(z^2 -1)*skew 
774
   h = h + (1/24)*(z^3 - 3*z)*exkurt - (1/36)*(2*z^3 - 5*z)*skew^2; 
775
   I1 = Ipower(1,h); I2 = Ipower(2,h); I3 = Ipower(3,h); I4 = Ipower(4,h);  I6 = Ipower(6,h); 
776

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

779
   E = dnorm(h)
780
   E = E + (1/24)*(   I4 - 6*I2 + 3*dnorm(h)   )*exkurt
781
   E = E +  (1/6)*(   I3 - 3*I1   )*skew;
782
   E = E + (1/72)*( I6 -15*I4+ 45*I2 - 15*dnorm(h) )*(skew^2)
783
   E = E/alpha
784

785
   MES = - location - sqrt(pm2)*min(-E,h)
786

787
   if(-E<=h){
788
      derMES = -mu + 0.5*(dpm2/sqrt(pm2))*E
789
      derE = (1/24)*(   I4 - 6*I2 + 3*dnorm(h)   )*derexkurt 
790
      derE = derE +  (1/6)*(   I3 - 3*I1   )*derskew
791
      derE = derE + (1/36)*(  I6 -15*I4 + 45*I2 - 15*dnorm(h) )*skew*derskew
792
      X = -h*dnorm(h) + (1/24)*(  derIpower(4,h) - 6*derIpower(2,h) -3*h*dnorm(h)  )*exkurt 
793
      X = X + (1/6)*( derIpower(3,h) - 3*derIpower(1,h) )*skew 
794
      X = X + (1/72)*( derIpower(6,h) - 15*derIpower(4,h) + 45*derIpower(2,h) + 15*h*dnorm(h)  )*skew^2
795
      derE = derE+derh*X  # X is a scalar
796
      derE = derE/alpha
797
      derMES = derMES + sqrt(pm2)*derE
798
    }else{
799
       derMES = -mu - 0.5*(dpm2/sqrt(pm2))*h - sqrt(pm2)*derh ; 
800
    }
801
    return( derMES )
802
}
803

804
centeredmoment = function(series,power)
805
{
806
   location = mean(series);
807
   out = sum( (series-location)^power  )/length(series);
808
   return(out);
809
}
810

811
operMES =  function(series,alpha,r)
812
{
813
   z = qnorm(alpha)
814
   location = mean(series);
815
   m2 = centeredmoment(series,2)
816
   m3 = centeredmoment(series,3)
817
   m4 = centeredmoment(series,4)
818
   skew = m3 / m2^(3/2);
819
   exkurt = m4 / m2^(2) - 3;
820

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

824
   MES = dnorm(h)
825
   MES = MES + (1/24)*(   Ipower(4,h) - 6*Ipower(2,h) + 3*dnorm(h)   )*exkurt
826
   MES = MES +  (1/6)*(   Ipower(3,h) - 3*Ipower(1,h)   )*skew;
827
   MES = MES + (1/72)*(  Ipower(6,h) -15*Ipower(4,h)+ 45*Ipower(2,h) - 15*dnorm(h) )*(skew^2)
828
   MES = - location - (sqrt(m2))*min( -MES/alpha , h )
829
   return(MES)
830
}
831

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

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

839
   # noincs - no of divisions in the axis labels
840
   # marks - type of marker for each y
841
   # legpos - legend position
842
   # leglabs - legend labels
843
  
844
   # rescale to fit on same axis
845
   scaledy2var <- (y2var - min(y2var)) / (max(y2var) - min(y2var))
846
   scaledy2var <- (scaledy2var * (max(y1var) - min(y1var))) + min(y1var)
847

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

852
   # make up some labels and positions
853
   y1labs <- round(seq(min(y1var), max(y1var), length=noincs),2)
854

855
   # convert these to the y2 axis scaling
856
   y2labs <- (y1labs - min(y1var)) / (max(y1var) - min(y1var))
857
   y2labs <- (y2labs * (max(y2var) - min(y2var))) + min(y2var)
858
   y2labs <- round(y2labs, 2)
859

860
   axis(1)
861
   axis(2, at=y1labs, labels=y1labs)
862
   axis(4, at=y1labs, labels=y2labs)
863
   mtext(labels[3], side=4, line=2)
864
   mtext(labels[2], side=2, line=2)
865
   box()
866

867
   legend( legend=leglabs, lty = c(1,3), bty="o", x=legpos)
868
}
869

870

871

872
Product

Resources

Company
