CoCalc -- estimators.R

GitHub Repository: braverock/portfolioanalytics
Path: blob/master/sandbox/paper_analysis/insample/estimators.R
¹⁴³³ views
1
# Function for data cleaning
2
#-------------------------------------------
3

4
cleaning = function(mData, alpha , trim=1e-3)
5
{
6
   T=dim(mData)[1];
7
   N=dim(mData)[2];
8
   MCD = covMcd(mData,alpha=1-alpha)
9
   mu = as.matrix(MCD$raw.center) #no reweighting
10
   sigma = MCD$raw.cov 
11
   invSigma = solve(sigma);
12
   vd2t = c();
13
   cleaneddata = mData
14
   outlierdate = c()
15

16
   # 1. Sort the data in function of their extremeness
17
   # Extremeness is proxied by the robustly estimated squared Mahalanbobis distance
18

19
   for(t in c(1:T) )
20
   {
21
      d2t = as.matrix(mData[t,]-mu)%*%invSigma%*%t(as.matrix(mData[t,]-mu));      
22
      vd2t = c(vd2t,d2t);
23
   }
24
   sortt = sort(vd2t,index.return=TRUE)$ix
25

26
   # 2. Outlier detection
27
   # 2.1. Only the alpha most extreme observations can be qualified as outliers
28

29
   T.alpha = floor(T * (1-alpha))+1
30
   cleanedt=sortt[c(T.alpha:T)]  
31

32
   # 2.2. multivariate winsorization (Khan et al, 2007) : 
33
   # bound the MD of the most exteme observations to a quantile of the chi squared distribution, N df
34

35
   for(t in cleanedt ){
36
        if(vd2t[t]>qchisq(1-trim,N)){
37
              # print(c("Observation",as.character(date[t]),"is detected as outlier and cleaned") );
38
               cleaneddata[t,] = sqrt(qchisq(1-trim,N)/vd2t[t])*mData[t,]; 
39
               outlierdate = c(outlierdate,date[t]) } }
40
   return(list(cleaneddata,outlierdate))
41
}
42

43
# Definition of statistics needed to compute Gaussian and modified VaR and ES for univariate time series
44
#---------------------------------------------------------------------------------------------------------
45

46
centeredmoment = function(series,power)
47
{
48
   location = mean(series);
49
   out = sum( (series-location)^power  )/length(series);
50
   return(out);
51
}
52

53
skew = function(series)
54
{
55
   m2 = centeredmoment(series,2)
56
   m3 = centeredmoment(series,3)
57
   out = m3 / m2^(3/2)
58
   return(out)
59
}
60

61
exkur =  function(series)
62
{
63
   m2 = centeredmoment(series,2)
64
   m4 = centeredmoment(series,4)
65
   out = m4 / m2^(2) - 3
66
   return(out)
67
}
68

69
pvalJB = function(series)
70
{
71
   m2 = centeredmoment(series,2)
72
   m3 = centeredmoment(series,3)
73
   m4 = centeredmoment(series,4)
74
   skew = m3 / m2^(3/2);
75
   exkur = m4 / m2^(2) - 3;
76
   JB = ( length(series)/6 )*( skew^2 + (1/4)*(exkur^2) )
77
   out = 1-pchisq(JB,df=2)
78
}
79

80
gausVaR =  function(series,alpha)
81
{
82
   location = mean(series);
83
   m2 = centeredmoment(series,2)
84
   out = - location - qnorm(alpha)*sqrt(m2)
85
   return(out)
86
}
87

88
gausES =  function(series,alpha)
89
{
90
   location = mean(series);
91
   m2 = centeredmoment(series,2)
92
   out = - location + dnorm(qnorm(alpha))*sqrt(m2)/alpha
93
   return(out)
94
}
95

96
MVaR =  function(series,alpha,r)
97
{
98
   z = qnorm(alpha)
99
   location = mean(series);
100
   m2 = centeredmoment(series,2)
101
   m3 = centeredmoment(series,3)
102
   m4 = centeredmoment(series,4)
103
   skew = m3 / m2^(3/2);
104
   exkurt = m4 / m2^(2) - 3;
105

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

109
   out = - location - h*sqrt(m2)
110
   return(out)
111
}
112

113
Ipower = function(power,h)
114
{
115

116
   fullprod = 1;
117

118
   if( (power%%2)==0 ) #even number: number mod is zero
119
   {
120
      pstar = power/2;
121
      for(j in c(1:pstar))
122
      {
123
         fullprod = fullprod*(2*j)
124
      }
125
      I = fullprod*dnorm(h);
126

127
      for(i in c(1:pstar) )
128
      { 
129
         prod = 1;
130
         for(j in c(1:i) )
131
         {
132
            prod = prod*(2*j)
133
         }
134
         I = I + (fullprod/prod)*(h^(2*i))*dnorm(h)
135
      }
136
   }
137
   else{
138
      pstar = (power-1)/2
139
      for(j in c(0:pstar) )
140
      {
141
         fullprod = fullprod*( (2*j)+1 )
142
      }
143
      I = -fullprod*pnorm(h);
144

145
      for(i in c(0:pstar) )
146
      { 
147
         prod = 1;
148
         for(j in c(0:i) )
149
         {
150
            prod = prod*( (2*j) + 1 )
151
         }
152
         I = I + (fullprod/prod)*(h^(  (2*i) + 1))*dnorm(h)
153
      }
154
   }
155
   return(I)
156
}
157

158
MES =  function(series,alpha,r=2)
159
{
160
   z = qnorm(alpha)
161
   location = mean(series);
162
   m2 = centeredmoment(series,2)
163
   m3 = centeredmoment(series,3)
164
   m4 = centeredmoment(series,4)
165
   skew = m3 / m2^(3/2);
166
   exkurt = m4 / m2^(2) - 3;
167

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

171
   MES = dnorm(h)
172
   MES = MES + (1/24)*(   Ipower(4,h) - 6*Ipower(2,h) + 3*dnorm(h)   )*exkurt
173
   MES = MES +  (1/6)*(   Ipower(3,h) - 3*Ipower(1,h)   )*skew;
174
   MES = MES + (1/72)*(  Ipower(6,h) -15*Ipower(4,h)+ 45*Ipower(2,h) - 15*dnorm(h) )*(skew^2)
175
   MES = - location + (sqrt(m2)/alpha)*MES
176
   return(MES)
177
}
178

179
operMES =  function(series,alpha,r)
180
{
181
   z = qnorm(alpha)
182
   location = mean(series);
183
   m2 = centeredmoment(series,2)
184
   m3 = centeredmoment(series,3)
185
   m4 = centeredmoment(series,4)
186
   skew = m3 / m2^(3/2);
187
   exkurt = m4 / m2^(2) - 3;
188

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

192
   MES = dnorm(h)
193
   MES = MES + (1/24)*(   Ipower(4,h) - 6*Ipower(2,h) + 3*dnorm(h)   )*exkurt
194
   MES = MES +  (1/6)*(   Ipower(3,h) - 3*Ipower(1,h)   )*skew;
195
   MES = MES + (1/72)*(  Ipower(6,h) -15*Ipower(4,h)+ 45*Ipower(2,h) - 15*dnorm(h) )*(skew^2)
196
   MES = - location - (sqrt(m2))*min( -MES/alpha , h )
197
   return(MES)
198
}
199

200
# Definition of statistics needed to compute Gaussian and modified VaR and ES for the return series of portfolios
201
# and to compute the contributions to portfolio downside risk, made by the different positions in the portfolio.
202
#----------------------------------------------------------------------------------------------------------------
203

204

205

206
m2 = function(w,sigma)
207
{
208
   return(t(w)%*%sigma%*%w)
209
}
210
derm2 = function(w,sigma)
211
{
212
   return(2*sigma%*%w)
213
}
214
m3 = function(w,M3)
215
{
216
   return(t(w)%*%M3%*%(w%x%w))
217
}
218
derm3 = function(w,M3)
219
{
220
   return(3*M3%*%(w%x%w))
221
}
222
m4 = function(w,M4)
223
{
224
   return(t(w)%*%M4%*%(w%x%w%x%w))
225
}
226
derm4 = function(w,M4)
227
{
228
   return(4*M4%*%(w%x%w%x%w))
229
}
230

231
precision = 4;
232

233
Portmean = function(w,mu,precision=4)
234
{
235
   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) )
236
}
237

238
Portsd =  function(w,sigma,precision=4)
239
{
240
   pm2 = m2(w,sigma)
241
   dpm2 = derm2(w,sigma)
242
   dersd = (0.5*as.vector(dpm2))/sqrt(pm2);
243
   contrib = dersd*as.vector(w)
244
   return(list(  round( sqrt(pm2) , precision ) , round( contrib , precision ) , round ( contrib/sqrt(pm2) , precision) )) 
245
}
246

247
PortgausVaR =  function(alpha,w,mu,sigma,precision=4)
248
{
249
   location = t(w)%*%mu
250
   pm2 = m2(w,sigma)
251
   dpm2 = derm2(w,sigma)
252
   VaR = - location - qnorm(alpha)*sqrt(pm2)
253
   derVaR = - as.vector(mu)- qnorm(alpha)*(0.5*as.vector(dpm2))/sqrt(pm2);
254
   contrib = derVaR*as.vector(w) 
255
   return(list( round( VaR , precision ) , round ( contrib , precision ) , round( contrib/VaR , precision) )) 
256
}
257

258
PortgausES =  function(alpha,w,mu,sigma,precision=4)
259
{
260
   location = t(w)%*%mu
261
   pm2 = m2(w,sigma)
262
   dpm2 = derm2(w,sigma)
263
   ES = - location + dnorm(qnorm(alpha))*sqrt(pm2)/alpha
264
   derES = - as.vector(mu) + (1/alpha)*dnorm(qnorm(alpha))*(0.5*as.vector(dpm2))/sqrt(pm2);
265
   contrib = as.vector(w)*derES;
266
   return(list( round( ES , precision ) , round( contrib , precision) , round( contrib/ES , precision) )) 
267
}
268

269
PortSkew =  function(w,sigma,M3)
270
{
271
   pm2 = m2(w,sigma)
272
   pm3 = m3(w,M3)
273
   skew = pm3 / pm2^(3/2);
274
   return( skew )
275
}
276

277
PortKurt =  function(w,sigma,M4)
278
{
279
   pm2 = m2(w,sigma)
280
   pm4 = m4(w,M4)
281
   kurt = pm4 / pm2^(2) ;
282
   return( kurt )
283
}
284

285
PortMVaR =  function(alpha,w,mu,sigma,M3,M4,precision=4)
286
{
287
   z = qnorm(alpha)
288
   location = t(w)%*%mu
289
   pm2 = m2(w,sigma)
290
   dpm2 = as.vector( derm2(w,sigma) )
291
   pm3 = m3(w,M3)
292
   dpm3 = as.vector( derm3(w,M3) )
293
   pm4 = m4(w,M4)
294
   dpm4 = as.vector( derm4(w,M4) )
295

296
   skew = pm3 / pm2^(3/2);
297
   exkurt = pm4 / pm2^(2) - 3;
298

299
   derskew = ( 2*(pm2^(3/2))*dpm3 - 3*pm3*sqrt(pm2)*dpm2 )/(2*pm2^3)
300
   derexkurt = ( (pm2)*dpm4 - 2*pm4*dpm2    )/(pm2^3)
301

302
   h = z + (1/6)*(z^2 -1)*skew 
303
   h = h + (1/24)*(z^3 - 3*z)*exkurt - (1/36)*(2*z^3 - 5*z)*skew^2; 
304

305
   MVaR = - location - h*sqrt(pm2)
306

307
   derGausVaR = - as.vector(mu)- qnorm(alpha)*(0.5*as.vector(dpm2))/sqrt(pm2);
308
   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 )
309
   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  )
310
   contrib = as.vector(w)*as.vector(derMVaR)
311
   return(list(  round( MVaR , precision) , round( contrib , precision ), round (contrib/MVaR , precision ) ) ) 
312
}
313

314
resPortMVaR =  function(alpha,w,mu,sigma,resSigma,M3,M4,precision=4)
315
{
316
   z = qnorm(alpha)
317
   location = t(w)%*%mu
318
   pm2 = m2(w,sigma)
319
   respm2 = m2(w,resSigma)
320
   dpm2 = as.vector( derm2(w,sigma) )
321
   drespm2 = as.vector( derm2(w,resSigma) )
322
   pm3 = m3(w,M3)
323
   dpm3 = as.vector( derm3(w,M3) )
324
   pm4 = m4(w,M4)
325
   dpm4 = as.vector( derm4(w,M4) )
326

327
   skew = pm3 / respm2^(3/2);
328
   exkurt = pm4 / respm2^(2) - 3;
329

330
   derskew = ( 2*(respm2^(3/2))*dpm3 - 3*pm3*sqrt(respm2)*drespm2 )/(2*respm2^3)
331
   derexkurt = ( (respm2)*dpm4 - 2*pm4*drespm2    )/(respm2^3)
332

333
   h = z + (1/6)*(z^2 -1)*skew 
334
   h = h + (1/24)*(z^3 - 3*z)*exkurt - (1/36)*(2*z^3 - 5*z)*skew^2; 
335

336
   MVaR = - location - h*sqrt(pm2)
337

338
   derGausVaR = - as.vector(mu)- qnorm(alpha)*(0.5*as.vector(dpm2))/sqrt(pm2);
339
   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 )
340
   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  )
341
   contrib = as.vector(w)*as.vector(derMVaR)
342
   return(list(  round( MVaR , precision) , round( contrib , precision ), round (contrib/MVaR , precision ) ) ) 
343
}
344

345
derIpower = function(power,h)
346
{
347

348
   fullprod = 1;
349

350
   if( (power%%2)==0 ) #even number: number mod is zero
351
   {
352
      pstar = power/2;
353
      for(j in c(1:pstar))
354
      {
355
         fullprod = fullprod*(2*j)
356
      }
357
      I = -fullprod*h*dnorm(h);
358

359
      for(i in c(1:pstar) )
360
      { 
361
         prod = 1;
362
         for(j in c(1:i) )
363
         {
364
            prod = prod*(2*j)
365
         }
366
         I = I + (fullprod/prod)*(h^(2*i-1))*(2*i-h^2)*dnorm(h)
367
      }
368
   }else{
369
      pstar = (power-1)/2
370
      for(j in c(0:pstar) )
371
      {
372
         fullprod = fullprod*( (2*j)+1 )
373
      }
374
      I = -fullprod*dnorm(h);
375

376
      for(i in c(0:pstar) )
377
      { 
378
         prod = 1;
379
         for(j in c(0:i) )
380
         {
381
            prod = prod*( (2*j) + 1 )
382
         }
383
         I = I + (fullprod/prod)*(h^(2*i)*(2*i+1-h^2) )*dnorm(h)
384
      }
385
   }
386
   return(I)
387
}
388

389

390
PortMES =  function(alpha,w,mu,sigma,M3,M4,precision=4)
391
{
392
   z = qnorm(alpha)
393
   location = t(w)%*%mu
394
   pm2 = m2(w,sigma)
395
   dpm2 = as.vector( derm2(w,sigma) )
396
   pm3 = m3(w,M3)
397
   dpm3 = as.vector( derm3(w,M3) )
398
   pm4 = m4(w,M4)
399
   dpm4 = as.vector( derm4(w,M4) )
400

401
   skew = pm3 / pm2^(3/2);
402
   exkurt = pm4 / pm2^(2) - 3;
403

404
   derskew = ( 2*(pm2^(3/2))*dpm3 - 3*pm3*sqrt(pm2)*dpm2 )/(2*pm2^3)
405
   derexkurt = ( (pm2)*dpm4 - 2*pm4*dpm2    )/(pm2^3)
406

407
   h = z + (1/6)*(z^2 -1)*skew 
408
   h = h + (1/24)*(z^3 - 3*z)*exkurt - (1/36)*(2*z^3 - 5*z)*skew^2; 
409

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

412
   E = dnorm(h)
413
   E = E + (1/24)*(   Ipower(4,h) - 6*Ipower(2,h) + 3*dnorm(h)   )*exkurt
414
   E = E +  (1/6)*(   Ipower(3,h) - 3*Ipower(1,h)   )*skew;
415
   E = E + (1/72)*(  Ipower(6,h) -15*Ipower(4,h)+ 45*Ipower(2,h) - 15*dnorm(h) )*(skew^2)
416
   E = E/alpha
417
   MES = - location + sqrt(pm2)*E
418

419
   derMES = -mu + 0.5*(dpm2/sqrt(pm2))*E
420
   derE = (1/24)*(   Ipower(4,h) - 6*Ipower(2,h) + 3*dnorm(h)   )*derexkurt 
421
   derE = derE +  (1/6)*(   Ipower(3,h) - 3*Ipower(1,h)   )*derskew
422
   derE = derE + (1/36)*(  Ipower(6,h) -15*Ipower(4,h)+ 45*Ipower(2,h) - 15*dnorm(h) )*skew*derskew
423
   X = -h*dnorm(h) + (1/24)*(  derIpower(4,h) - 6*derIpower(2,h) -3*h*dnorm(h)  )*exkurt 
424
   X = X + (1/6)*( derIpower(3,h) - 3*derIpower(1,h) )*skew 
425
   X = X + (1/72)*( derIpower(6,h) - 15*derIpower(4,h) + 45*derIpower(2,h) + 15*h*dnorm(h)  )*skew^2
426
   derE = derE+derh*X  # X is a scalar
427
   derE = derE/alpha
428
   derMES = derMES + sqrt(pm2)*derE
429
   contrib = as.vector(w)*as.vector(derMES)
430
   return(list(  round( MES , precision ) , round( contrib , precision ), round( contrib/MES, precision )) ) 
431
}
432

433

434
operPortMES =  function(alpha,w,mu,sigma,M3,M4,precision=4)
435
{
436
   z = qnorm(alpha)
437
   location = t(w)%*%mu
438
   pm2 = m2(w,sigma)
439
   dpm2 = as.vector( derm2(w,sigma) )
440
   pm3 = m3(w,M3)
441
   dpm3 = as.vector( derm3(w,M3) )
442
   pm4 = m4(w,M4)
443
   dpm4 = as.vector( derm4(w,M4) )
444

445
   skew = pm3 / pm2^(3/2);
446
   exkurt = pm4 / pm2^(2) - 3;
447

448
   derskew = ( 2*(pm2^(3/2))*dpm3 - 3*pm3*sqrt(pm2)*dpm2 )/(2*pm2^3)
449
   derexkurt = ( (pm2)*dpm4 - 2*pm4*dpm2    )/(pm2^3)
450

451
   h = z + (1/6)*(z^2 -1)*skew 
452
   h = h + (1/24)*(z^3 - 3*z)*exkurt - (1/36)*(2*z^3 - 5*z)*skew^2; 
453
   I1 = Ipower(1,h); I2 = Ipower(2,h); I3 = Ipower(3,h); I4 = Ipower(4,h);  I6 = Ipower(6,h); 
454

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

457
   E = dnorm(h)
458
   E = E + (1/24)*(   I4 - 6*I2 + 3*dnorm(h)   )*exkurt
459
   E = E +  (1/6)*(   I3 - 3*I1   )*skew;
460
   E = E + (1/72)*( I6 -15*I4+ 45*I2 - 15*dnorm(h) )*(skew^2)
461
   E = E/alpha
462

463
   MES = - location - sqrt(pm2)*min(-E,h)
464

465
   if(-E<=h){
466
      derMES = -mu + 0.5*(dpm2/sqrt(pm2))*E
467
      derE = (1/24)*(   I4 - 6*I2 + 3*dnorm(h)   )*derexkurt 
468
      derE = derE +  (1/6)*(   I3 - 3*I1   )*derskew
469
      derE = derE + (1/36)*(  I6 -15*I4 + 45*I2 - 15*dnorm(h) )*skew*derskew
470
      X = -h*dnorm(h) + (1/24)*(  derIpower(4,h) - 6*derIpower(2,h) -3*h*dnorm(h)  )*exkurt 
471
      X = X + (1/6)*( derIpower(3,h) - 3*derIpower(1,h) )*skew 
472
      X = X + (1/72)*( derIpower(6,h) - 15*derIpower(4,h) + 45*derIpower(2,h) + 15*h*dnorm(h)  )*skew^2
473
      derE = derE+derh*X  # X is a scalar
474
      derE = derE/alpha
475
      derMES = derMES + sqrt(pm2)*derE }else{
476
      derMES = -mu - 0.5*(dpm2/sqrt(pm2))*h - sqrt(pm2)*derh ;  }
477
   contrib = as.vector(w)*as.vector(derMES)
478
   return(list( round( MES, precision) , round( contrib , precision ) , round(contrib/MES,precision) ) ) 
479
}
480

481

482
resoperPortMES =  function(alpha,w,mu,sigma,resSigma,M3,M4,precision=4)
483
{
484
   z = qnorm(alpha)
485
   location = t(w)%*%mu
486
   pm2 = m2(w,sigma)
487
   respm2 = m2(w,resSigma)
488
   dpm2 = as.vector( derm2(w,sigma) )
489
   drespm2 = as.vector( derm2(w,resSigma) )
490
   pm3 = m3(w,M3)
491
   dpm3 = as.vector( derm3(w,M3) )
492
   pm4 = m4(w,M4)
493
   dpm4 = as.vector( derm4(w,M4) )
494

495
   skew = pm3 / respm2^(3/2);
496
   exkurt = pm4 / respm2^(2) - 3;
497

498
   derskew = ( 2*(respm2^(3/2))*dpm3 - 3*pm3*sqrt(respm2)*drespm2 )/(2*respm2^3)
499
   derexkurt = ( (respm2)*dpm4 - 2*pm4*drespm2    )/(respm2^3)
500

501
   h = z + (1/6)*(z^2 -1)*skew 
502
   h = h + (1/24)*(z^3 - 3*z)*exkurt - (1/36)*(2*z^3 - 5*z)*skew^2; 
503
   I1 = Ipower(1,h); I2 = Ipower(2,h); I3 = Ipower(3,h); I4 = Ipower(4,h);  I6 = Ipower(6,h); 
504

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

507
   E = dnorm(h)
508
   E = E + (1/24)*(   I4 - 6*I2 + 3*dnorm(h)   )*exkurt
509
   E = E +  (1/6)*(   I3 - 3*I1   )*skew;
510
   E = E + (1/72)*( I6 -15*I4+ 45*I2 - 15*dnorm(h) )*(skew^2)
511
   E = E/alpha
512

513
   MES = - location - sqrt(pm2)*min(-E,h)
514

515
   if(-E<=h){
516
      derMES = -mu + 0.5*(dpm2/sqrt(pm2))*E
517
      derE = (1/24)*(   I4 - 6*I2 + 3*dnorm(h)   )*derexkurt 
518
      derE = derE +  (1/6)*(   I3 - 3*I1   )*derskew
519
      derE = derE + (1/36)*(  I6 -15*I4 + 45*I2 - 15*dnorm(h) )*skew*derskew
520
      X = -h*dnorm(h) + (1/24)*(  derIpower(4,h) - 6*derIpower(2,h) -3*h*dnorm(h)  )*exkurt 
521
      X = X + (1/6)*( derIpower(3,h) - 3*derIpower(1,h) )*skew 
522
      X = X + (1/72)*( derIpower(6,h) - 15*derIpower(4,h) + 45*derIpower(2,h) + 15*h*dnorm(h)  )*skew^2
523
      derE = derE+derh*X  # X is a scalar
524
      derE = derE/alpha
525
      derMES = derMES + sqrt(pm2)*derE }else{
526
      derMES = -mu - 0.5*(dpm2/sqrt(pm2))*h - sqrt(pm2)*derh ;  }
527
   contrib = as.vector(w)*as.vector(derMES)
528
   return(list( round( MES, precision) , round( contrib , precision ) , round(contrib/MES,precision) ) ) 
529
}
530

531
realportES = function(data,w,VaR,precision=4)
532
{
533
    T = dim(data)[1]
534
    N = dim(data)[2]
535
    c_exceed = 0;
536
    r_exceed = 0;
537
    realizedcontrib = rep(0,N);
538
    for(t in c(1:T) )
539
    {
540
       rt = as.vector(data[t,])
541
       rp = sum(w*rt)
542
       if(rp<=-VaR){
543
          c_exceed = c_exceed + 1;
544
          r_exceed = r_exceed + rp;
545
          for( i in c(1:N) ){
546
             realizedcontrib[i] =realizedcontrib[i] + w[i]*rt[i] }
547
       }
548
    }
549
    realizedcontrib=as.numeric(realizedcontrib)/r_exceed ; 
550
    return( list(round(-r_exceed/c_exceed,precision),c_exceed,round(realizedcontrib,precision)) )
551
}
552

553
realportVaR = function(data,w,alpha,precision=4)
554
{
555
    portret = c();
556
    T = dim(data)[1]
557
    N = dim(data)[2]
558
    for( t in c(1:T) ){
559
       portret = c(portret,sum(w*as.numeric(data[t,])))
560
    }
561
    VaR = sort(portret)[floor(alpha*T)]
562
    return(-VaR)
563
}
564
Product

Resources

Company
