1
2 #include <R.h>
3 #include <Rinternals.h>
4
5 // // //
6 // // Sample estimators
7 // // //
8
M3sample(SEXP XX,SEXP NN,SEXP PP,SEXP CC)9 SEXP M3sample(SEXP XX, SEXP NN, SEXP PP, SEXP CC){
10 /*
11 arguments
12 XX : vector of length NN * PP with observations
13 NN : integer, number of observations
14 PP : integer, number of assets
15 CC : numeric, normalising constant
16
17 Written by Dries Cornilly
18 */
19
20 // // declare pointers for the vector
21 double *X;
22
23 // use REAL() to access the C array inside the numeric vector passed in from R
24 X = REAL(XX);
25 int N = asInteger(NN);
26 int P = asInteger(PP);
27 double c = asReal(CC);
28
29 // allocate and compute the coskewness matrix
30 SEXP M3S = PROTECT(allocVector(REALSXP, P * (P + 1) * (P + 2) / 6));
31 double *rM3S = REAL(M3S);
32
33 int iter = 0;
34 // loop over the unique indices i <= j <= k
35 for (int ii = 0; ii < P; ii++) {
36 int iiN = ii * N;
37 for (int jj = ii; jj < P; jj++) {
38 int jjN = jj * N;
39 for (int kk = jj; kk < P; kk++) {
40 int kkN = kk * N;
41 // compute coskewness element
42 double elem = 0.0;
43 for (int tt = 0; tt < N; tt++) {
44 elem += X[iiN + tt] * X[jjN + tt] * X[kkN + tt];
45 }
46
47 // fill in the element and update the iterator
48 rM3S[iter] = c * elem;
49 iter++;
50 } // loop kk
51 } // loop jj
52 } // loop ii
53
54 UNPROTECT(1);
55 return M3S;
56 }
57
M4sample(SEXP XX,SEXP NN,SEXP PP)58 SEXP M4sample(SEXP XX, SEXP NN, SEXP PP){
59 /*
60 arguments
61 XX : vector of length NN * PP with observations
62 NN : integer, number of observations
63 PP : integer, number of assets
64
65 Written by Dries Cornilly
66 */
67
68 // // declare pointers for the vectors
69 double *X;
70
71 // use REAL() to access the C array inside the numeric vector passed in from R
72 X = REAL(XX);
73 int N = asInteger(NN);
74 int P = asInteger(PP);
75 double n = asReal(NN);
76
77 // allocate and compute the cokurtosis matrix
78 SEXP M4S = PROTECT(allocVector(REALSXP, P * (P + 1) * (P + 2) * (P + 3) / 24));
79 double *rM4S = REAL(M4S);
80
81 int iter = 0;
82 // compute the unique elements
83 for (int ii = 0; ii < P; ii++) {
84 int iiN = ii * N;
85 for (int jj = ii; jj < P; jj++) {
86 int jjN = jj * N;
87 for (int kk = jj; kk < P; kk++) {
88 int kkN = kk * N;
89 for (int ll = kk; ll < P; ll++) {
90 int llN = ll * N;
91 // compute cokurtosis element
92 double elem = 0.0;
93 for (int tt = 0; tt < N; tt++) {
94 elem += X[iiN + tt] * X[jjN + tt] * X[kkN + tt] * X[llN + tt];
95 }
96
97 // fill in the element and update the iterator
98 rM4S[iter] = elem / n;
99 iter++;
100 } // loop ll
101 } // loop kk
102 } // loop jj
103 } // loop ii
104
105 UNPROTECT(1);
106 return M4S;
107 }
108
109
110 // // //
111 // // Shrinkage and structured estimators for covariance
112 // // //
113
VM2(SEXP mm11,SEXP mm22,SEXP NN,SEXP PP)114 SEXP VM2(SEXP mm11, SEXP mm22, SEXP NN, SEXP PP){
115 /*
116 arguments
117 mm11 : numeric vector of t(Xc) %*% Xc / NN
118 mm22 : numeric vector of t(Xc^2) %*% Xc^2 / NN
119 NN : integer, number of observations
120 PP : integer, number of assets
121
122 Written by Dries Cornilly
123 */
124
125 // // declare pointers for the vector
126 double *m11, *m22;
127
128 // use REAL() to access the C array inside the numeric vector passed in from R
129 m11 = REAL(mm11);
130 m22 = REAL(mm22);
131 int P = asInteger(PP);
132 double N = asReal(NN);
133
134 // allocate and compute the sample variance and covariances
135 SEXP VM2vec = PROTECT(allocVector(REALSXP, 3));
136 double *rVM2vec = REAL(VM2vec);
137 for (int ii = 0; ii < 3; ii++) {
138 rVM2vec[ii] = 0.0;
139 }
140
141 // loop over the unique indices i <= j
142 for (int ii = 0; ii < P; ii++) {
143 int iiP = ii * P;
144 for (int jj = ii; jj < P; jj++) {
145 int jjP = jj * P;
146 if (ii == jj) {
147 // element s_ii
148 double temp = m22[iiP + ii] - m11[iiP + ii] * m11[iiP + ii];
149 rVM2vec[0] += temp / N;
150 rVM2vec[2] += temp / N;
151 } else {
152 // element s_ij
153 rVM2vec[0] += 2.0 * (m22[jjP + ii] - m11[jjP + ii] * m11[jjP + ii]) / N;
154 }
155 } // loop jj
156 } // loop ii
157
158 // compute cov(T1, Sigma)
159 rVM2vec[1] = rVM2vec[2];
160
161 for (int ii = 0; ii < P; ii++) {
162 int iiP = ii * P;
163 for (int jj = ii + 1; jj < P; jj++) {
164 int jjP = jj * P;
165 rVM2vec[1] += 2.0 * (m22[jjP + ii] - m11[iiP + ii] * m11[jjP + jj]) / N;
166 } // loop over jj
167 } // loop over ii
168 rVM2vec[1] /= P;
169
170 UNPROTECT(1);
171 return VM2vec;
172 }
173
CM2_1F(SEXP XXc,SEXP ffcobs,SEXP ffvar,SEXP mm11,SEXP mm22,SEXP NN,SEXP PP)174 SEXP CM2_1F(SEXP XXc, SEXP ffcobs, SEXP ffvar, SEXP mm11, SEXP mm22, SEXP NN, SEXP PP){
175 /*
176 arguments
177 XXc : numeric vector with centered observations
178 ffcobs : numeric vector with centered factor observations
179 ffvar : double with factor variance
180 mm11 : numeric vector of t(Xc) %*% Xc / NN
181 mm22 : numeric vector of t(Xc^2) %*% Xc^2 / NN
182 NN : integer, number of observations
183 PP : integer, number of assets
184
185 Written by Dries Cornilly
186 */
187
188 // // declare pointers for the vector
189 double *Xc, *m11, *m22, *fcobs;
190
191 // use REAL() to access the C array inside the numeric vector passed in from R
192 Xc = REAL(XXc);
193 m11 = REAL(mm11);
194 m22 = REAL(mm22);
195 fcobs = REAL(ffcobs);
196 double fvar = asReal(ffvar);
197 double N = asReal(NN);
198 int n = asInteger(NN);
199 int P = asInteger(PP);
200
201 // compute some helper variables for later on
202 double fvar2 = fvar * fvar;
203 double covXf[P];
204 for (int ii = 0; ii < P; ii++) {
205 int iiN = ii * n;
206 covXf[ii] = 0.0;
207 for (int tt = 0; tt < n; tt++) {
208 covXf[ii] += Xc[iiN + tt] * fcobs[tt];
209 }
210 covXf[ii] /= N;
211 }
212
213 // allocate and compute the sample variance and covariances
214 SEXP CM2vec = PROTECT(allocVector(REALSXP, 1));
215 double *rCM2vec = REAL(CM2vec);
216 rCM2vec[0] = 0.0;
217
218 // loop over the unique indices i <= j
219 for (int ii = 0; ii < P; ii++) {
220 int iiP = ii * P;
221 int iiN = ii * n;
222 for (int jj = ii; jj < P; jj++) {
223 int jjP = jj * P;
224 int jjN = jj * n;
225 if (ii == jj) {
226 // element s_ii
227 rCM2vec[0] += (m22[iiP + ii] - m11[iiP + ii] * m11[iiP + ii]) / N;
228 } else {
229 // element s_ij
230 double S211 = 0.0;
231 double S121 = 0.0;
232 double S112 = 0.0;
233 for (int tt = 0; tt < n; tt++) {
234 S211 += Xc[iiN + tt] * Xc[iiN + tt] * Xc[jjN + tt] * fcobs[tt];
235 S121 += Xc[iiN + tt] * Xc[jjN + tt] * Xc[jjN + tt] * fcobs[tt];
236 S112 += Xc[iiN + tt] * Xc[jjN + tt] * fcobs[tt] * fcobs[tt];
237 }
238
239 double temp_i0 = S211 / N - m11[jjP + ii] * covXf[ii];
240 double temp_j0 = S121 / N - m11[jjP + ii] * covXf[jj];
241 double temp_var = S112 / N - m11[jjP + ii] * fvar;
242
243 rCM2vec[0] += 2.0 * (covXf[jj] * temp_i0 / fvar + covXf[ii] * temp_j0 / fvar -
244 covXf[ii] * covXf[jj] * temp_var / fvar2) / N;
245 }
246 } // loop jj
247 } // loop ii
248
249 UNPROTECT(1);
250 return CM2vec;
251 }
252
CM2_CC(SEXP XXc,SEXP rrcoef,SEXP mm11,SEXP mm22,SEXP NN,SEXP PP)253 SEXP CM2_CC(SEXP XXc, SEXP rrcoef, SEXP mm11, SEXP mm22, SEXP NN, SEXP PP){
254 /*
255 arguments
256 XXc : numeric vector with centered observations
257 rrcoef : double with factor variance
258 mm11 : numeric vector of t(Xc) %*% Xc / NN
259 mm22 : numeric vector of t(Xc^2) %*% Xc^2 / NN
260 NN : integer, number of observations
261 PP : integer, number of assets
262
263 Written by Dries Cornilly
264 */
265
266 // // declare pointers for the vector
267 double *Xc, *m11, *m22;
268
269 // use REAL() to access the C array inside the numeric vector passed in from R
270 Xc = REAL(XXc);
271 m11 = REAL(mm11);
272 m22 = REAL(mm22);
273 double rcoef = asReal(rrcoef);
274 double N = asReal(NN);
275 int n = asInteger(NN);
276 int P = asInteger(PP);
277
278 // allocate and compute the sample variance and covariances
279 SEXP CM2vec = PROTECT(allocVector(REALSXP, 1));
280 double *rCM2vec = REAL(CM2vec);
281 rCM2vec[0] = 0.0;
282
283 // loop over the unique indices i <= j
284 for (int ii = 0; ii < P; ii++) {
285 int iiP = ii * P;
286 int iiN = ii * n;
287 for (int jj = ii; jj < P; jj++) {
288 int jjP = jj * P;
289 int jjN = jj * n;
290 if (ii == jj) {
291 // element s_ii
292 rCM2vec[0] += (m22[iiP + ii] - m11[iiP + ii] * m11[iiP + ii]) / N;
293 } else {
294 // element s_ij
295 double S31 = 0.0;
296 double S13 = 0.0;
297 for (int tt = 0; tt < n; tt++) {
298 S31 += Xc[iiN + tt] * Xc[iiN + tt] * Xc[iiN + tt] * Xc[jjN + tt];
299 S13 += Xc[iiN + tt] * Xc[jjN + tt] * Xc[jjN + tt] * Xc[jjN + tt];
300 }
301
302 double temp_ii = S31 / N - m11[jjP + ii] * m11[iiP + ii];
303 double temp_jj = S13 / N - m11[jjP + ii] * m11[jjP + jj];
304
305 rCM2vec[0] += rcoef * (sqrt(m11[jjP + jj] / m11[iiP + ii]) * temp_ii +
306 sqrt(m11[iiP + ii] / m11[jjP + jj]) * temp_jj) / N;
307 }
308 } // loop jj
309 } // loop ii
310
311 UNPROTECT(1);
312 return CM2vec;
313 }
314
315
316 // // //
317 // // Shrinkage and structured estimators for coskewness
318 // // //
319
M3_T23(SEXP mmargskews,SEXP PP)320 SEXP M3_T23(SEXP mmargskews, SEXP PP){
321 /*
322 arguments
323 mmargskews : vector of length PP with marginal skewness values
324 PP : integer, number of assets
325
326 Written by Dries Cornilly
327 */
328
329 // // declare pointers for the vector
330 double *margskews;
331
332 // use REAL() to access the C array inside the numeric vector passed in from R
333 margskews = REAL(mmargskews);
334 int P = asInteger(PP);
335
336 // allocate and compute the coskewness matrix
337 SEXP M3 = PROTECT(allocVector(REALSXP, P * (P + 1) * (P + 2) / 6));
338 double *rM3 = REAL(M3);
339
340 int iter = 0;
341 // loop over the unique indices i <= j <= k
342 for (int ii = 0; ii < P; ii++) {
343 for (int jj = ii; jj < P; jj++) {
344 for (int kk = jj; kk < P; kk++) {
345 // compute coskewness element
346 double elem = 0.0;
347 if ((ii == jj) && (jj == kk)) elem = margskews[ii];
348
349 // fill in the element and update the iterator
350 rM3[iter] = elem;
351 iter++;
352 } // loop kk
353 } // loop jj
354 } // loop ii
355
356 UNPROTECT(1);
357 return M3;
358 }
359
M3_Simaan(SEXP mmargskewsroot,SEXP PP)360 SEXP M3_Simaan(SEXP mmargskewsroot, SEXP PP){
361 /*
362 arguments
363 mmargskewsroot : vector of length PP with marginal skewness values (phi^(1 / 3))
364 PP : integer, number of assets
365
366 Written by Dries Cornilly, coskewness matrix based on Simaan (1993)
367 */
368
369 // // declare pointers for the vector
370 double *margskewsroot;
371
372 // use REAL() to access the C array inside the numeric vector passed in from R
373 margskewsroot = REAL(mmargskewsroot);
374 int P = asInteger(PP);
375
376 // allocate and compute the coskewness matrix
377 SEXP M3 = PROTECT(allocVector(REALSXP, P * (P + 1) * (P + 2) / 6));
378 double *rM3 = REAL(M3);
379
380 int iter = 0;
381 // loop over the unique indices i <= j <= k
382 for (int ii = 0; ii < P; ii++) {
383 for (int jj = ii; jj < P; jj++) {
384 for (int kk = jj; kk < P; kk++) {
385 // compute and fill coskewness element
386 rM3[iter] = margskewsroot[ii] * margskewsroot[jj] * margskewsroot[kk];
387 iter++;
388 } // loop kk
389 } // loop jj
390 } // loop ii
391
392 UNPROTECT(1);
393 return M3;
394 }
395
M3_1F(SEXP mmargskews,SEXP bbeta,SEXP ffskew,SEXP PP)396 SEXP M3_1F(SEXP mmargskews, SEXP bbeta, SEXP ffskew, SEXP PP){
397 /*
398 arguments
399 mmargskews : vector of length PP with marginal skewness values
400 bbeta : vector of length PP with factor loadings
401 ffskew : factor skewness
402 PP : integer, number of assets
403
404 Written by Dries Cornilly
405 */
406
407 // // declare pointers for the vector
408 double *margskews, *beta;
409
410 // use REAL() to access the C array inside the numeric vector passed in from R
411 margskews = REAL(mmargskews);
412 beta = REAL(bbeta);
413 double fskew = asReal(ffskew);
414 int P = asInteger(PP);
415
416 // allocate and compute the coskewness matrix
417 SEXP M3 = PROTECT(allocVector(REALSXP, P * (P + 1) * (P + 2) / 6));
418 double *rM3 = REAL(M3);
419
420 int iter = 0;
421 // loop over the unique indices i <= j <= k
422 for (int ii = 0; ii < P; ii++) {
423 for (int jj = ii; jj < P; jj++) {
424 for (int kk = jj; kk < P; kk++) {
425 // compute and fill coskewness element
426 if ((ii == jj) && (jj == kk)) {
427 rM3[iter] = margskews[ii];
428 } else {
429 rM3[iter] = beta[ii] * beta[jj] * beta[kk] * fskew;
430 }
431 iter++;
432 } // loop kk
433 } // loop jj
434 } // loop ii
435
436 UNPROTECT(1);
437 return M3;
438 }
439
M3_CCoefficients(SEXP mmargvars,SEXP mmargkurts,SEXP mm21,SEXP mm22,SEXP XXc,SEXP NN,SEXP PP)440 SEXP M3_CCoefficients(SEXP mmargvars, SEXP mmargkurts,
441 SEXP mm21, SEXP mm22, SEXP XXc, SEXP NN, SEXP PP){
442 /*
443 arguments
444 mmargvars : vector of length PP with marginal variances
445 mmargkurts : vector of length PP with marginal kurtosis values
446 mm21 : t(Xc^2) %*% Xc
447 mm22 : t(Xc^2) %*% Xc^2
448 XXc : numeric vector with centered observations
449 NN : integer, number of observations
450 PP : integer, number of assets
451
452 Written by Dries Cornilly
453 */
454
455 // // declare pointers for the vector
456 double *margvars, *margkurts, *m21, *m22, *Xc;
457
458 // use REAL() to access the C array inside the numeric vector passed in from R
459 margvars = REAL(mmargvars);
460 margkurts = REAL(mmargkurts);
461 m21 = REAL(mm21);
462 m22 = REAL(mm22);
463 Xc = REAL(XXc);
464 double N = asReal(NN);
465 int n = asInteger(NN);
466 int P = asInteger(PP);
467 double p = asReal(PP);
468
469 // allocate and compute the extended correlation coefficients
470 SEXP CCoef = PROTECT(allocVector(REALSXP, 3));
471 double *rCCoef = REAL(CCoef);
472
473 // compute the generalized correlation coefficients
474 rCCoef[0] = 0.0;
475 rCCoef[2] = 0.0;
476 for (int ii = 0; ii < P; ii++) {
477 for (int jj = ii + 1; jj < P; jj++) {
478 int jjP = jj * P;
479 rCCoef[0] += m21[jjP + ii] / sqrt(margkurts[ii] * margvars[jj]);
480 rCCoef[2] += m22[jjP + ii] / sqrt(margkurts[ii] * margkurts[jj]);
481 }
482 }
483 rCCoef[0] *= 2.0 / (p * (p - 1.0));
484 rCCoef[2] *= 2.0 / (p * (p - 1.0));
485 rCCoef[1] = 0.0;
486 for (int ii = 0; ii < P; ii++) {
487 int iiN = ii * n;
488 for (int jj = ii + 1; jj < P; jj++) {
489 int jjN = jj * n;
490 for (int kk = jj + 1; kk < P; kk++) {
491 int kkN = kk * n;
492 double m111 = 0.0;
493 for (int tt = 0; tt < n; tt++) {
494 m111 += Xc[iiN + tt] * Xc[jjN + tt] * Xc[kkN + tt];
495 }
496 m111 /= N;
497 double nc = sqrt(margvars[ii] * rCCoef[2] * sqrt(margkurts[jj] * margkurts[kk])) +
498 sqrt(margvars[jj] * rCCoef[2] * sqrt(margkurts[ii] * margkurts[kk])) +
499 sqrt(margvars[kk] * rCCoef[2] * sqrt(margkurts[ii] * margkurts[jj]));
500 rCCoef[1] += m111 / nc;
501 }
502 }
503 }
504 rCCoef[1] *= 6.0 / (p * (p - 1.0) * (p - 2.0));
505
506 UNPROTECT(1);
507 return CCoef;
508 }
509
M3_CC(SEXP mmargvars,SEXP mmargskews,SEXP mmargkurts,SEXP rr2,SEXP rr4,SEXP rr5,SEXP PP)510 SEXP M3_CC(SEXP mmargvars, SEXP mmargskews, SEXP mmargkurts,
511 SEXP rr2, SEXP rr4, SEXP rr5, SEXP PP){
512 /*
513 arguments
514 mmargvars : vector of length PP with marginal variances
515 mmargskews : vector of length PP with marginal skewness values
516 mmargkurts : vector of length PP with marginal kurtosis values
517 rr2 : generalized correlation of order 2
518 rr4 : generalized correlation of order 4
519 rr5 : generalized correlation of order 5
520 PP : integer, number of assets
521
522 Written by Dries Cornilly
523 */
524
525 // // declare pointers for the vector
526 double *margvars, *margskews, *margkurts;
527
528 // use REAL() to access the C array inside the numeric vector passed in from R
529 margvars = REAL(mmargvars);
530 margskews = REAL(mmargskews);
531 margkurts = REAL(mmargkurts);
532 double r2 = asReal(rr2);
533 double r4 = asReal(rr4);
534 double r5 = asReal(rr5);
535 int P = asInteger(PP);
536
537 // allocate and compute the coskewness matrix
538 SEXP M3 = PROTECT(allocVector(REALSXP, P * (P + 1) * (P + 2) / 6));
539 double *rM3 = REAL(M3);
540
541 int iter = 0;
542 // loop over the unique indices i <= j <= k
543 for (int ii = 0; ii < P; ii++) {
544 for (int jj = ii; jj < P; jj++) {
545 for (int kk = jj; kk < P; kk++) {
546 // compute and fill coskewness element
547 if (ii == jj) {
548 if (jj == kk) {
549 // phi_iii
550 rM3[iter] = margskews[ii];
551 } else {
552 // phi_iik
553 rM3[iter] = r2 * sqrt(margkurts[kk] * margvars[ii]);
554 }
555 } else {
556 if (jj == kk) {
557 // phi_ijj
558 rM3[iter] = r2 * sqrt(margvars[ii] * margkurts[jj]);
559 } else {
560 // phi_ijk
561 rM3[iter] = r4 * sqrt(r5) * (sqrt(margvars[kk] * sqrt(margkurts[ii] * margkurts[jj])) +
562 sqrt(margvars[jj] * sqrt(margkurts[ii] * margkurts[kk])) +
563 sqrt(margvars[ii] * sqrt(margkurts[jj] * margkurts[kk]))) / 3.0;
564 }
565 }
566 iter++;
567 } // loop kk
568 } // loop jj
569 } // loop ii
570
571 UNPROTECT(1);
572 return M3;
573 }
574
VM3kstat(SEXP XXc,SEXP XXc2,SEXP SS11,SEXP SS21,SEXP SS22,SEXP SS31,SEXP SS42,SEXP SS33,SEXP NN,SEXP PP)575 SEXP VM3kstat(SEXP XXc, SEXP XXc2,
576 SEXP SS11, SEXP SS21, SEXP SS22, SEXP SS31,
577 SEXP SS42, SEXP SS33, SEXP NN, SEXP PP){
578 /*
579 arguments
580 XXc : numeric vector with centered observations
581 XXc2 : numeric vector with centered and squared observations
582 SS11 : numeric vector of t(Xc) %*% Xc
583 SS21 : numeric vector of t(Xc^2) %*% Xc
584 SS22 : numeric vector of t(Xc^2) %*% Xc^2
585 SS31 : numeric vector of t(Xc^3) %*% Xc
586 SS42 : numeric vector of t(Xc^4) %*% Xc^2
587 SS33 : numeric vector of t(Xc^3) %*% Xc^3
588 NN : integer, number of observations
589 PP : integer, number of assets
590
591 Written by Dries Cornilly
592 */
593
594 // // declare pointers for the vector
595 double *Xc, *Xc2, *S11, *S21, *S22, *S31, *S42, *S33;
596
597 // use REAL() to access the C array inside the numeric vector passed in from R
598 Xc = REAL(XXc);
599 Xc2 = REAL(XXc2);
600 S11 = REAL(SS11);
601 S21 = REAL(SS21);
602 S22 = REAL(SS22);
603 S31 = REAL(SS31);
604 S42 = REAL(SS42);
605 S33 = REAL(SS33);
606 double N = asReal(NN);
607 int n = asInteger(NN);
608 int P = asInteger(PP);
609
610 // initialize useful constants
611 double N2 = N * N;
612 double N3 = N2 * N;
613 double N4 = N3 * N;
614 double N5 = N4 * N;
615 double N122 = (N - 1.0) * (N - 1.0) * (N - 2.0) * (N - 2.0);
616
617 double alpha = N * N122 * (N - 3.0) * (N - 4.0) * (N - 5.0);
618 double ciii_S2_3 = (9.0 * N4 - 72.0 * N3 + 213.0 * N2 - 270.0 * N + 120.0) / alpha;
619 double ciii_S4S2 = (-6.0 * N5 + 33.0 * N4 - 42.0 * N3 - 75.0 * N2 + 210.0 * N - 120.0) / alpha;
620 double ciii_S3_2 = (-N5 - 4.0 * N4 + 41.0 * N3 - 40.0 * N2 - 100.0 * N + 80.0) / alpha;
621 double ciii_S6 = (N5 - 5.0 * N4 + 13.0 * N3 - 23.0 * N2 + 22.0 * N - 8.0) /
622 (N122 * (N - 3.0) * (N - 4.0) * (N - 5.0));
623
624 double ciij_S02S20_2 = (N4 - 8.0 * N3 + 29.0 * N2 - 46.0 * N + 24.0) / alpha;
625 double ciij_S20S11_2 = (8.0 * N4 - 64.0 * N3 + 184.0 * N2 - 224.0 * N + 96.0) / alpha;
626 double ciij_S20S22 = (-2.0 * N5 + 10.0 * N4 - 10.0 * N3 - 34.0 * N2 + 84.0 * N - 48.0) / alpha;
627 double ciij_S31S11 = (-4.0 * N5 + 24.0 * N4 - 36.0 * N3 - 32.0 * N2 + 112.0 * N - 64.0) / alpha;
628 double ciij_S40S02 = (-N4 + 4.0 * N3 - 9.0 * N2 + 14.0 * N - 8.0) / alpha;
629 double ciij_S21_2 = (-N5 + 25.0 * N3 - 36.0 * N2 - 60.0 * N + 48.0) / alpha;
630 double ciij_S12S30 = (-4.0 * N4 + 16.0 * N3 - 4.0 * N2 - 40.0 * N + 32.0) / alpha;
631
632 double cijk_S002S110_2 = (N4 - 8.0 * N3 + 25.0 * N2 - 34.0 * N + 16.0) / alpha;
633 double cijk_S011S101S110 = (6.0 * N4 - 48.0 * N3 + 134.0 * N2 - 156.0 * N + 64.0) / alpha;
634 double cijk_S112S110 = (-2.0 * N5 + 12.0 * N4 - 18.0 * N3 - 16.0 * N2 + 56.0 * N - 32.0) / alpha;
635 double cijk_S002S020S200 = (4.0 * N2 - 12.0 * N + 8.0) / alpha;
636 double cijk_S022S200 = (-N4 + 4.0 * N3 - 9.0 * N2 + 14.0 * N - 8.0) / alpha;
637 double cijk_S111_2 = (-N5 + 2.0 * N4 + 17.0 * N3 - 34.0 * N2 - 40.0 * N + 32.0) / alpha;
638 double cijk_S102S120 = (-2.0 * N4 + 8.0 * N3 - 2.0 * N2 - 20.0 * N + 16.0) / alpha;
639
640 // allocate and compute the sample variance and covariances
641 SEXP VM3vec = PROTECT(allocVector(REALSXP, 3));
642 double *rVM3vec = REAL(VM3vec);
643 for (int ii = 0; ii < 3; ii++) {
644 rVM3vec[ii] = 0.0;
645 }
646
647 // loop over the unique indices i <= j <= k
648 for (int ii = 0; ii < P; ii++) {
649 int iiP = ii * P;
650 int iiN = ii * n;
651 for (int jj = ii; jj < P; jj++) {
652 int jjP = jj * P;
653 int jjN = jj * n;
654 for (int kk = jj; kk < P; kk++) {
655 int kkP = kk * P;
656 int kkN = kk * n;
657 if (ii == jj) {
658 if (jj == kk) {
659 // element phi_iii
660 double temp = ciii_S2_3 * S11[iiP + ii] * S11[iiP + ii] * S11[iiP + ii] +
661 ciii_S4S2 * S22[iiP + ii]* S11[iiP + ii] +
662 ciii_S3_2 * S21[iiP + ii] * S21[iiP + ii] + ciii_S6 * S42[iiP + ii];
663 rVM3vec[0] += temp;
664 rVM3vec[2] += temp;
665 } else {
666 // element phi_iik
667 rVM3vec[0] += 3.0 * (ciij_S02S20_2 * S11[kkP + kk] * S11[iiP + ii] * S11[iiP + ii] +
668 ciij_S20S11_2 * S11[iiP + ii] * S11[kkP + ii] * S11[kkP + ii] +
669 ciij_S20S22 * S11[iiP + ii] * S22[kkP + ii] + ciij_S31S11 * S31[kkP + ii] * S11[kkP + ii] +
670 ciij_S40S02 * S22[iiP + ii] * S11[kkP + kk] + ciij_S21_2 * S21[kkP + ii] * S21[kkP + ii] +
671 ciij_S12S30 * S21[iiP + kk] * S21[iiP + ii] + ciii_S6 * S42[kkP + ii]);
672 }
673 } else {
674 if (jj == kk) {
675 // element phi_ijj
676 rVM3vec[0] += 3.0 * (ciij_S02S20_2 * S11[iiP + ii] * S11[jjP + jj] * S11[jjP + jj] +
677 ciij_S20S11_2 * S11[jjP + jj] * S11[iiP + jj] * S11[iiP + jj] +
678 ciij_S20S22 * S11[jjP + jj] * S22[iiP + jj] + ciij_S31S11 * S31[iiP + jj] * S11[iiP + jj] +
679 ciij_S40S02 * S22[jjP + jj] * S11[iiP + ii] + ciij_S21_2 * S21[iiP + jj] * S21[iiP + jj] +
680 ciij_S12S30 * S21[jjP + ii] * S21[jjP + jj] + ciii_S6 * S42[iiP + jj]);
681 } else {
682 // element phi_ijk
683 double S211i = 0.0;
684 double S211j = 0.0;
685 double S211k = 0.0;
686 double S111 = 0.0;
687 double S222 = 0.0;
688
689 for (int tt = 0; tt < n; tt++) {
690 S211i += Xc2[iiN + tt] * Xc[jjN + tt] * Xc[kkN + tt];
691 S211j += Xc[iiN + tt] * Xc2[jjN + tt] * Xc[kkN + tt];
692 S211k += Xc[iiN + tt] * Xc[jjN + tt] * Xc2[kkN + tt];
693 S111 += Xc[iiN + tt] * Xc[jjN + tt] * Xc[kkN + tt];
694 S222 += Xc2[iiN + tt] * Xc2[jjN + tt] * Xc2[kkN + tt];
695 }
696
697 rVM3vec[0] += 6.0 * (cijk_S002S110_2 * S11[kkP + kk] * S11[jjP + ii] * S11[jjP + ii] +
698 (cijk_S011S101S110 * S11[kkP + jj] * S11[kkP + ii] + cijk_S112S110 * S211k) * S11[jjP + ii] +
699 cijk_S002S110_2 * S11[jjP + jj] * S11[kkP + ii] * S11[kkP + ii] +
700 cijk_S112S110 * S211j * S11[kkP + ii] +
701 cijk_S002S110_2 * S11[iiP + ii] * S11[kkP + jj] * S11[kkP + jj] +
702 cijk_S112S110 * S211i * S11[kkP + jj] +
703 (cijk_S002S020S200 * S11[jjP + jj] * S11[kkP + kk] + cijk_S022S200 * S22[kkP + jj]) * S11[iiP + ii] +
704 cijk_S022S200 * S22[kkP + ii] * S11[jjP + jj] + cijk_S022S200 * S22[jjP + ii] * S11[kkP + kk] +
705 cijk_S111_2 * S111 * S111 +
706 cijk_S102S120 * S21[iiP + kk] * S21[iiP + jj] + cijk_S102S120 * S21[jjP + kk] * S21[jjP + ii] +
707 cijk_S102S120 * S21[kkP + jj] * S21[kkP + ii] + ciii_S6 * S222);
708 }
709 }
710 } // loop kk
711 } // loop jj
712 } // loop ii
713
714 // compute cov(T1, Phi)
715 rVM3vec[1] = rVM3vec[2];
716
717 double c_S11_3 = (24.0 * N2 - 72.0 * N + 48.0) / alpha;
718 double c_S02S20S11 = (9.0 * N4 - 72.0 * N3 + 189.0 * N2 - 198.0 * N + 72.0) / alpha;
719 double c_S22S11 = (-9.0 * N4 + 36.0 * N3 - 81.0 * N2 + 126.0 * N - 72.0) / alpha;
720 double c_S13S20 = (-3.0 * N5 + 21.0 * N4 - 39.0 * N3 + 3.0 * N2 + 42.0 * N - 24.0) / alpha;
721 double c_S21S12 = (-9.0 * N4 + 36.0 * N3 - 9.0 * N2 - 90.0 * N + 72.0) / alpha;
722 double c_S03S30 = (-N5 + 5.0 * N4 + 5.0 * N3 - 31.0 * N2 - 10.0 * N + 8.0) / alpha;
723
724 for (int ii = 0; ii < P; ii++) {
725 int iiP = ii * P;
726 for (int jj = ii + 1; jj < P; jj++) {
727 int jjP = jj * P;
728 rVM3vec[1] += 2 * (c_S11_3 * S11[jjP + ii] * S11[jjP + ii] * S11[jjP + ii] +
729 (c_S02S20S11 * S11[iiP + ii] * S11[jjP + jj] + c_S22S11 * S22[jjP + ii]) * S11[jjP + ii] +
730 c_S13S20 * S31[iiP + jj] * S11[iiP + ii] + c_S13S20 * S31[jjP + ii] * S11[jjP + jj] +
731 c_S21S12 * S21[jjP + ii] * S21[iiP + jj] + c_S03S30 * S21[iiP + ii] * S21[jjP + jj] +
732 ciii_S6 * S33[jjP + ii]);
733 } // loop over jj
734 } // loop over ii
735 rVM3vec[1] /= P;
736
737 UNPROTECT(1);
738 return VM3vec;
739 }
740
VM3(SEXP XXc,SEXP XXc2,SEXP mm11,SEXP mm21,SEXP mm22,SEXP mm31,SEXP mm42,SEXP mm33,SEXP NN,SEXP PP)741 SEXP VM3(SEXP XXc, SEXP XXc2, SEXP mm11, SEXP mm21, SEXP mm22, SEXP mm31,
742 SEXP mm42, SEXP mm33, SEXP NN, SEXP PP){
743 /*
744 arguments
745 XXc : numeric vector with centered observations
746 XXc2 : numeric vector with centered and squared observations
747 mm11 : numeric vector of t(Xc) %*% Xc / NN
748 mm21 : numeric vector of t(Xc^2) %*% Xc / NN
749 mm22 : numeric vector of t(Xc^2) %*% Xc^2 / NN
750 mm31 : numeric vector of t(Xc^3) %*% Xc / NN
751 mm42 : numeric vector of t(Xc^4) %*% Xc^2 / NN
752 mm33 : numeric vector of t(Xc^3) %*% Xc^3 / NN
753 NN : integer, number of observations
754 PP : integer, number of assets
755
756 Written by Dries Cornilly
757 */
758
759 // // declare pointers for the vector
760 double *Xc, *Xc2, *m11, *m21, *m22, *m31, *m42, *m33;
761
762 // use REAL() to access the C array inside the numeric vector passed in from R
763 Xc = REAL(XXc);
764 Xc2 = REAL(XXc2);
765 m11 = REAL(mm11);
766 m21 = REAL(mm21);
767 m22 = REAL(mm22);
768 m31 = REAL(mm31);
769 m42 = REAL(mm42);
770 m33 = REAL(mm33);
771 double N = asReal(NN);
772 int n = asInteger(NN);
773 int P = asInteger(PP);
774
775 // allocate and compute the sample variance and covariances
776 SEXP VM3vec = PROTECT(allocVector(REALSXP, 3));
777 double *rVM3vec = REAL(VM3vec);
778 for (int ii = 0; ii < 3; ii++) {
779 rVM3vec[ii] = 0.0;
780 }
781
782 // loop over the unique indices i <= j <= k
783 for (int ii = 0; ii < P; ii++) {
784 int iiP = ii * P;
785 int iiN = ii * n;
786 for (int jj = ii; jj < P; jj++) {
787 int jjP = jj * P;
788 int jjN = jj * n;
789 for (int kk = jj; kk < P; kk++) {
790 int kkP = kk * P;
791 int kkN = kk * n;
792 if (ii == jj) {
793 if (jj == kk) {
794 // element phi_iii
795 double temp = (m42[iiP + ii] - m21[iiP + ii] * m21[iiP + ii] - 6.0 * m22[iiP + ii] * m11[iiP + ii] +
796 9.0 * m11[iiP + ii] * m11[iiP + ii] * m11[iiP + ii]) / N;
797 rVM3vec[0] += temp;
798 rVM3vec[2] += temp;
799 } else {
800 // element phi_iik
801 rVM3vec[0] += 3.0 * (m42[kkP + ii] - m21[kkP + ii] * m21[kkP + ii] - 4.0 * m31[kkP + ii] * m11[kkP + ii] -
802 2.0 * m22[kkP + ii] * m11[iiP + ii] + 8.0 * m11[iiP + ii] * m11[kkP + ii] * m11[kkP + ii] +
803 m11[kkP + kk] * m11[iiP + ii] * m11[iiP + ii]) / N;
804 }
805 } else {
806 if (jj == kk) {
807 // element phi_ijj
808 rVM3vec[0] += 3.0 * (m42[iiP + jj] - m21[iiP + jj] * m21[iiP + jj] - 4.0 * m31[iiP + jj] * m11[iiP + jj] -
809 2.0 * m22[iiP + jj] * m11[jjP + jj] + 8.0 * m11[jjP + jj] * m11[iiP + jj] * m11[iiP + jj] +
810 m11[iiP + ii] * m11[jjP + jj] * m11[jjP + jj]) / N;
811 } else {
812 // element phi_ijk
813 double S211 = 0.0;
814 double S121 = 0.0;
815 double S112 = 0.0;
816 double S111 = 0.0;
817 double S222 = 0.0;
818
819 for (int tt = 0; tt < n; tt++) {
820 S211 += Xc2[iiN + tt] * Xc[jjN + tt] * Xc[kkN + tt];
821 S211 += Xc[iiN + tt] * Xc2[jjN + tt] * Xc[kkN + tt];
822 S211 += Xc[iiN + tt] * Xc[jjN + tt] * Xc2[kkN + tt];
823 S111 += Xc[iiN + tt] * Xc[jjN + tt] * Xc[kkN + tt];
824 S222 += Xc2[iiN + tt] * Xc2[jjN + tt] * Xc2[kkN + tt];
825 }
826
827 rVM3vec[0] += 6.0 * (S222 / N - S111 * S111 / (N * N) - 2.0 * S211 / N * m11[kkP + jj] -
828 2.0 * S121 / N * m11[kkP + ii] - 2.0 * S112 / N * m11[jjP + ii] +
829 6.0 * m11[kkP + ii] * m11[kkP + jj] * m11[jjP + ii] +
830 m11[iiP + ii] * m11[kkP + jj] * m11[kkP + jj] + m11[jjP + jj] * m11[kkP + ii] * m11[kkP + ii] +
831 m11[kkP + kk] * m11[jjP + ii] * m11[jjP + ii]) / N;
832 }
833 }
834 } // loop kk
835 } // loop jj
836 } // loop ii
837
838 // compute cov(T1, Phi)
839 rVM3vec[1] = rVM3vec[2];
840
841 for (int ii = 0; ii < P; ii++) {
842 int iiP = ii * P;
843 for (int jj = ii + 1; jj < P; jj++) {
844 int jjP = jj * P;
845 rVM3vec[1] += 2 * (m33[jjP + ii] - m21[iiP + ii] * m21[jjP + jj] - 3.0 * m31[jjP + ii] * m11[jjP + jj] -
846 3.0 * m31[iiP + jj] * m11[iiP + ii] + 9.0 * m11[iiP + ii] * m11[jjP + jj] * m11[jjP + ii]) / N;
847 } // loop over jj
848 } // loop over ii
849 rVM3vec[1] /= P;
850
851 UNPROTECT(1);
852 return VM3vec;
853 }
854
CM3_Simaan(SEXP XXc,SEXP XXc2,SEXP mmargskewsroot,SEXP mm11,SEXP mm21,SEXP mm22,SEXP mm31,SEXP mm42,SEXP mm51,SEXP NN,SEXP PP)855 SEXP CM3_Simaan(SEXP XXc, SEXP XXc2, SEXP mmargskewsroot,
856 SEXP mm11, SEXP mm21, SEXP mm22, SEXP mm31,
857 SEXP mm42, SEXP mm51, SEXP NN, SEXP PP){
858 /*
859 arguments
860 XXc : numeric vector with centered observations
861 XXc2 : numeric vector with centered and squared observations
862 mmargskewsroot : phi^(-2 / 3) with phi the marginal skewness values
863 mm11 : numeric vector of t(Xc) %*% Xc / NN
864 mm21 : numeric vector of t(Xc^2) %*% Xc / NN
865 mm22 : numeric vector of t(Xc^2) %*% Xc^2 / NN
866 mm31 : numeric vector of t(Xc^3) %*% Xc / NN
867 mm42 : numeric vector of t(Xc^4) %*% Xc^2 / NN
868 mm51 : numeric vector of t(Xc^5) %*% Xc / NN
869 NN : integer, number of observations
870 PP : integer, number of assets
871
872 Written by Dries Cornilly
873 */
874
875 // // declare pointers for the vector
876 double *Xc, *Xc2, *m11, *m21, *m22, *m31, *m42, *m51, *margskewsroot;
877
878 // use REAL() to access the C array inside the numeric vector passed in from R
879 Xc = REAL(XXc);
880 Xc2 = REAL(XXc2);
881 m11 = REAL(mm11);
882 m21 = REAL(mm21);
883 m22 = REAL(mm22);
884 m31 = REAL(mm31);
885 m42 = REAL(mm42);
886 m51 = REAL(mm51);
887 margskewsroot = REAL(mmargskewsroot);
888 double N = asReal(NN);
889 int n = asInteger(NN);
890 int P = asInteger(PP);
891
892 // allocate and compute the sample variance and covariances
893 SEXP CM3vec = PROTECT(allocVector(REALSXP, 1));
894 double *rCM3vec = REAL(CM3vec);
895 rCM3vec[0] = 0.0;
896
897 // loop over the unique indices i <= j <= k
898 for (int ii = 0; ii < P; ii++) {
899 int iiP = ii * P;
900 int iiN = ii * n;
901 for (int jj = ii; jj < P; jj++) {
902 int jjP = jj * P;
903 int jjN = jj * n;
904 for (int kk = jj; kk < P; kk++) {
905 int kkP = kk * P;
906 int kkN = kk * n;
907 if (ii == jj) {
908 if (jj == kk) {
909 // element phi_iii
910 rCM3vec[0] += (m42[iiP + ii] - m21[iiP + ii] * m21[iiP + ii] - 6.0 * m22[iiP + ii] * m11[iiP + ii] +
911 9.0 * m11[iiP + ii] * m11[iiP + ii] * m11[iiP + ii]);
912 } else {
913 // element phi_iik
914 double temp_ii = m51[kkP + ii] - m21[kkP + ii] * m21[iiP + ii] - 4.0 * m31[kkP + ii] * m11[iiP + ii] -
915 2.0 * m22[iiP + ii] * m11[kkP + ii] + 9.0 * m11[iiP + ii] * m11[iiP + ii] * m11[kkP + ii];
916 double temp_kk = m42[iiP + kk] - m21[kkP + ii] * m21[kkP + kk] - 3.0 * m22[kkP + ii] * m11[kkP + kk] -
917 m22[kkP + kk] * m11[iiP + ii] - 2.0 * m31[iiP + kk] * m11[kkP + ii] +
918 6.0 * m11[kkP + kk] * m11[kkP + ii] * m11[kkP + ii] + 3.0 * m11[kkP + kk] * m11[kkP + kk] * m11[iiP + ii];
919
920 rCM3vec[0] += margskewsroot[ii] * margskewsroot[ii] * margskewsroot[kk] *
921 (2.0 * m21[jjP + jj] * m21[kkP + kk] * temp_ii + m21[iiP + ii] * m21[iiP + ii] * temp_kk);
922 }
923 } else {
924 if (jj == kk) {
925 // element phi_ijj
926 double temp_ii = m42[jjP + ii] - m21[iiP + jj] * m21[iiP + ii] - 3.0 * m22[iiP + jj] * m11[iiP + ii] -
927 m22[iiP + ii] * m11[jjP + jj] - 2.0 * m31[jjP + ii] * m11[iiP + jj] +
928 6.0 * m11[iiP + ii] * m11[iiP + jj] * m11[iiP + jj] + 3.0 * m11[iiP + ii] * m11[iiP + ii] * m11[jjP + jj];
929 double temp_jj = m51[iiP + jj] - m21[iiP + jj] * m21[jjP + jj] - 4.0 * m31[iiP + jj] * m11[jjP + jj] -
930 2.0 * m22[jjP + jj] * m11[iiP + jj] + 9.0 * m11[jjP + jj] * m11[jjP + jj] * m11[iiP + jj];
931
932 rCM3vec[0] += margskewsroot[ii] * margskewsroot[jj] * margskewsroot[jj] *
933 (m21[jjP + jj] * m21[jjP + jj] * temp_ii + 2.0 * m21[iiP + ii] * m21[kkP + kk] * temp_jj);
934 } else {
935 // element phi_ijk
936 double S411 = 0.0;
937 double S141 = 0.0;
938 double S114 = 0.0;
939 double S111 = 0.0;
940 double S211 = 0.0;
941 double S121 = 0.0;
942 double S112 = 0.0;
943 for (int tt = 0; tt < n; tt++) {
944 S411 += Xc2[iiN + tt] * Xc2[iiN + tt] * Xc[jjN + tt] * Xc[kkN + tt];
945 S141 += Xc[iiN + tt] * Xc2[jjN + tt] * Xc2[jjN + tt] * Xc[kkN + tt];
946 S114 += Xc[iiN + tt] * Xc[jjN + tt] * Xc2[kkN + tt] * Xc2[kkN + tt];
947 S111 += Xc[iiN + tt] * Xc[jjN + tt] * Xc[kkN + tt];
948 S211 += Xc2[iiN + tt] * Xc[jjN + tt] * Xc[kkN + tt];
949 S121 += Xc[iiN + tt] * Xc2[jjN + tt] * Xc[kkN + tt];
950 S112 += Xc[iiN + tt] * Xc[jjN + tt] * Xc2[kkN + tt];
951 }
952
953 double temp_ii = S411 / N - S111 / N * m21[iiP + ii] - 3.0 * S211 / N * m11[iiP + ii] -
954 m22[iiP + ii] * m11[kkP + jj] - m31[jjP + ii] * m11[kkP + ii] - m31[kkP + ii] * m11[jjP + ii] +
955 6.0 * m11[iiP + ii] * m11[jjP + ii] * m11[kkP + ii] + 3.0 * m11[iiP + ii] * m11[iiP + ii] * m11[kkP + jj];
956 double temp_jj = S141 / N - S111 / N * m21[jjP + jj] - 3.0 * S121 / N * m11[jjP + jj] -
957 m22[jjP + jj] * m11[kkP + ii] - m31[iiP + jj] * m11[kkP + jj] - m31[kkP + jj] * m11[iiP + jj] +
958 6.0 * m11[jjP + jj] * m11[iiP + jj] * m11[kkP + jj] + 3.0 * m11[jjP + jj] * m11[jjP + jj] * m11[kkP + ii];
959 double temp_kk = S114 / N - S111 / N * m21[kkP + kk] - 3.0 * S112 / N * m11[kkP + kk] -
960 m22[kkP + kk] * m11[jjP + ii] - m31[jjP + kk] * m11[iiP + kk] - m31[iiP + kk] * m11[jjP + kk] +
961 6.0 * m11[kkP + kk] * m11[jjP + kk] * m11[iiP + kk] + 3.0 * m11[kkP + kk] * m11[kkP + kk] * m11[jjP + ii];
962
963 rCM3vec[0] += 2.0 * margskewsroot[ii] * margskewsroot[jj] * margskewsroot[kk] *
964 (m21[jjP + jj] * m21[kkP + kk] * temp_ii + m21[iiP + ii] * m21[kkP + kk] * temp_jj +
965 m21[iiP + ii] * m21[jjP + jj] * temp_kk);
966 }
967 }
968 } // loop kk
969 } // loop jj
970 } // loop ii
971
972 rCM3vec[0] /= N;
973
974 UNPROTECT(1);
975 return CM3vec;
976 }
977
CM3_1F(SEXP XXc,SEXP XXc2,SEXP ffcobs,SEXP ffvar,SEXP ffskew,SEXP mm11,SEXP mm21,SEXP mm22,SEXP mm42,SEXP NN,SEXP PP)978 SEXP CM3_1F(SEXP XXc, SEXP XXc2, SEXP ffcobs, SEXP ffvar, SEXP ffskew,
979 SEXP mm11, SEXP mm21, SEXP mm22, SEXP mm42, SEXP NN, SEXP PP){
980 /*
981 arguments
982 XXc : numeric vector with centered observations
983 XXc2 : numeric vector with centered and squared observations
984 ffcobs : numeric vector with centered factor observations
985 ffvar : double with the factor variance
986 ffskew : double with the factor skewness
987 mm11 : numeric vector of t(Xc) %*% Xc / NN
988 mm21 : numeric vector of t(Xc^2) %*% Xc / NN
989 mm22 : numeric vector of t(Xc^2) %*% Xc^2 / NN
990 mm42 : numeric vector of t(Xc^4) %*% Xc^2 / NN
991 NN : integer, number of observations
992 PP : integer, number of assets
993
994 Written by Dries Cornilly
995 */
996
997 // // declare pointers for the vector
998 double *Xc, *Xc2, *m11, *m21, *m22, *m42, *fcobs;
999
1000 // use REAL() to access the C array inside the numeric vector passed in from R
1001 Xc = REAL(XXc);
1002 Xc2 = REAL(XXc2);
1003 m11 = REAL(mm11);
1004 m21 = REAL(mm21);
1005 m22 = REAL(mm22);
1006 m42 = REAL(mm42);
1007 fcobs = REAL(ffcobs);
1008 double fvar = asReal(ffvar);
1009 double fskew = asReal(ffskew);
1010 double N = asReal(NN);
1011 int n = asInteger(NN);
1012 int P = asInteger(PP);
1013
1014 // compute some helper variables for later on
1015 double fvar3 = fvar * fvar * fvar;
1016 double covXf[P];
1017 double X1f2[P];
1018 double X1f3[P];
1019 double X11f1[P * P];
1020 for (int ii = 0; ii < P; ii++) {
1021 int iiN = ii * n;
1022 covXf[ii] = 0.0;
1023 X1f2[ii] = 0.0;
1024 X1f3[ii] = 0.0;
1025 for (int tt = 0; tt < n; tt++) {
1026 covXf[ii] += Xc[iiN + tt] * fcobs[tt];
1027 X1f2[ii] += Xc[iiN + tt] * fcobs[tt] * fcobs[tt];
1028 X1f3[ii] += Xc[iiN + tt] * fcobs[tt] * fcobs[tt] * fcobs[tt];
1029 }
1030 covXf[ii] /= N;
1031 X1f2[ii] /= N;
1032 X1f3[ii] /= N;
1033 for (int jj = ii; jj < P; jj++) {
1034 int jjN = jj * n;
1035 double elem = 0.0;
1036 for (int tt = 0; tt < n; tt++) {
1037 elem += Xc[iiN + tt] * Xc[jjN + tt] * fcobs[tt];
1038 }
1039 elem /= N;
1040 X11f1[ii * P + jj] = elem;
1041 X11f1[jj * P + ii] = elem;
1042 }
1043 }
1044
1045 // allocate and compute the sample variance and covariances
1046 SEXP CM3vec = PROTECT(allocVector(REALSXP, 1));
1047 double *rCM3vec = REAL(CM3vec);
1048 rCM3vec[0] = 0.0;
1049
1050 // loop over the unique indices i <= j <= k
1051 for (int ii = 0; ii < P; ii++) {
1052 int iiP = ii * P;
1053 int iiN = ii * n;
1054 for (int jj = ii; jj < P; jj++) {
1055 int jjP = jj * P;
1056 int jjN = jj * n;
1057 for (int kk = jj; kk < P; kk++) {
1058 int kkP = kk * P;
1059 int kkN = kk * n;
1060 if (ii == jj) {
1061 if (jj == kk) {
1062 // element phi_iii
1063 rCM3vec[0] += m42[iiP + ii] - m21[iiP + ii] * m21[iiP + ii] - 6.0 * m22[iiP + ii] * m11[iiP + ii] +
1064 9.0 * m11[iiP + ii] * m11[iiP + ii] * m11[iiP + ii];
1065 } else {
1066 // element phi_iik
1067 double S311 = 0.0;
1068 double S221 = 0.0;
1069 double S213 = 0.0;
1070 double S211 = 0.0;
1071 double S212 = 0.0;
1072 for (int tt = 0; tt < n; tt++) {
1073 S311 += Xc2[iiN + tt] * Xc[iiN + tt] * Xc[kkN + tt] * fcobs[tt];
1074 S221 += Xc2[iiN + tt] * Xc2[kkN + tt] * fcobs[tt];
1075 S213 += Xc2[iiN + tt] * Xc[kkN + tt] * fcobs[tt] * fcobs[tt] * fcobs[tt];
1076 S211 += Xc2[iiN + tt] * Xc[kkN + tt] * fcobs[tt];
1077 S212 += Xc2[iiN + tt] * Xc[kkN + tt] * fcobs[tt] * fcobs[tt];
1078 }
1079
1080 double temp_ii = S311 / N - m21[kkP + ii] * covXf[ii] - 2.0 * m11[kkP + ii] * X11f1[iiP + ii] -
1081 m11[iiP + ii] * X11f1[kkP + ii];
1082 double temp_kk = S221 / N - m21[kkP + ii] * covXf[kk] - m11[iiP + ii] * X11f1[kkP + kk] -
1083 2.0 * m11[kkP + ii] * X11f1[kkP + ii];
1084 double temp_fskew = S213 / N - m21[kkP + ii] * fskew - 3.0 * S211 / N * fvar -
1085 2.0 * X1f3[ii] * m11[kkP + ii] - X1f3[kk] * m11[iiP + ii] +
1086 3.0 * m11[iiP + ii] * covXf[kk] * fvar + 6.0 * m11[kkP + ii] * covXf[ii] * fvar;
1087 double temp_fvar = S212 / N - m21[kkP + ii] * fvar - 2.0 * X1f2[ii] * m11[kkP + ii] -
1088 X1f2[kk] * m11[iiP + ii];
1089
1090 rCM3vec[0] += 3.0 * ((2.0 * covXf[ii] * covXf[kk] * temp_ii +
1091 covXf[ii] * covXf[ii] * temp_kk) * fskew + covXf[ii] * covXf[ii] * covXf[kk] * temp_fskew -
1092 3.0 * covXf[ii] * covXf[ii] * covXf[kk] * fskew * temp_fvar / fvar) / fvar3;
1093 }
1094 } else {
1095 if (jj == kk) {
1096 // element phi_ijj
1097 double S221 = 0.0;
1098 double S131 = 0.0;
1099 double S123 = 0.0;
1100 double S121 = 0.0;
1101 double S122 = 0.0;
1102 for (int tt = 0; tt < n; tt++) {
1103 S221 += Xc2[iiN + tt] * Xc2[jjN + tt] * fcobs[tt];
1104 S131 += Xc[iiN + tt] * Xc2[jjN + tt] * Xc[jjN + tt] * fcobs[tt];
1105 S123 += Xc[iiN + tt] * Xc2[jjN + tt] * fcobs[tt] * fcobs[tt] * fcobs[tt];
1106 S121 += Xc[iiN + tt] * Xc2[jjN + tt] * fcobs[tt];
1107 S122 += Xc[iiN + tt] * Xc2[jjN + tt] * fcobs[tt] * fcobs[tt];
1108 }
1109 double temp_ii = S221 / N - m21[iiP + jj] * covXf[ii] - m11[jjP + jj] * X11f1[iiP + ii] -
1110 2.0 * m11[jjP + ii] * X11f1[jjP + ii];
1111 double temp_jj = S131 / N - m21[iiP + jj] * covXf[jj] - 2.0 * m11[jjP + ii] * X11f1[jjP + jj] -
1112 m11[jjP + jj] * X11f1[jjP + ii];
1113 double temp_fskew = S123 / N - m21[iiP + jj] * fskew - 3.0 * S121 / N * fvar -
1114 X1f3[ii] * m11[jjP + jj] - 2.0 * X1f3[jj] * m11[jjP + ii] +
1115 6.0 * m11[jjP + ii] * covXf[jj] * fvar + 3.0 * m11[jjP + jj] * covXf[ii] * fvar;
1116 double temp_fvar = S122 / N - m21[iiP + jj] * fvar - X1f2[ii] * m11[jjP + jj] -
1117 2.0 * X1f2[jj] * m11[jjP + ii];
1118
1119 rCM3vec[0] += 3.0 * ((covXf[jj] * covXf[jj] * temp_ii + 2.0 * covXf[ii] * covXf[jj] * temp_jj) * fskew +
1120 covXf[ii] * covXf[jj] * covXf[jj] * temp_fskew -
1121 3.0 * covXf[ii] * covXf[jj] * covXf[jj] * fskew * temp_fvar / fvar) / fvar3;
1122 } else {
1123 // element phi_ijk
1124 double S2111 = 0.0;
1125 double S1211 = 0.0;
1126 double S1121 = 0.0;
1127 double S1112 = 0.0;
1128 double S1111 = 0.0;
1129 double S1113 = 0.0;
1130 double S1110 = 0.0;
1131 for (int tt = 0; tt < n; tt++) {
1132 S2111 += Xc2[iiN + tt] * Xc[jjN + tt] * Xc[kkN + tt] * fcobs[tt];
1133 S1211 += Xc[iiN + tt] * Xc2[jjN + tt] * Xc[kkN + tt] * fcobs[tt];
1134 S1121 += Xc[iiN + tt] * Xc[jjN + tt] * Xc2[kkN + tt] * fcobs[tt];
1135 S1112 += Xc[iiN + tt] * Xc[jjN + tt] * Xc[kkN + tt] * fcobs[tt] * fcobs[tt];
1136 S1113 += Xc[iiN + tt] * Xc[jjN + tt] * Xc[kkN + tt] * fcobs[tt] * fcobs[tt] * fcobs[tt];
1137 S1111 += Xc[iiN + tt] * Xc[jjN + tt] * Xc[kkN + tt] * fcobs[tt];
1138 S1110 += Xc[iiN + tt] * Xc[jjN + tt] * Xc[kkN + tt];
1139 }
1140
1141 double temp_ii = S2111 / N - S1110 / N * covXf[ii] - m11[kkP + jj] * X11f1[iiP + ii] -
1142 m11[kkP + ii] * X11f1[jjP + ii] - m11[jjP + ii] * X11f1[kkP + ii];
1143 double temp_jj = S1211 / N - S1110 / N * covXf[jj] - m11[kkP + ii] * X11f1[jjP + jj] -
1144 m11[kkP + jj] * X11f1[iiP + jj] - m11[iiP + jj] * X11f1[kkP + jj];
1145 double temp_kk = S1121 / N - S1110 / N * covXf[kk] - m11[jjP + ii] * X11f1[kkP + kk] -
1146 m11[jjP + kk] * X11f1[iiP + kk] - m11[iiP + kk] * X11f1[jjP + kk];
1147 double temp_fskew = S1113 / N - S1110 / N * fskew - 3.0 * S1111 / N * fvar -
1148 X1f3[ii] * m11[kkP + jj] - X1f3[jj] * m11[kkP + ii] - X1f3[kk] * m11[jjP + ii] +
1149 3.0 * m11[jjP + ii] * covXf[kk] * fvar + 3.0 * m11[kkP + ii] * covXf[jj] * fvar +
1150 3.0 * m11[kkP + jj] * covXf[ii] * fvar;
1151 double temp_fvar = S1112 / N - S1110 / N * fvar - X1f2[ii] * m11[kkP + jj] -
1152 X1f2[jj] * m11[kkP + ii] - X1f2[kk] * m11[jjP + ii];
1153
1154 rCM3vec[0] += 6.0 * ((covXf[jj] * covXf[kk] * temp_ii + covXf[ii] * covXf[kk] * temp_jj +
1155 covXf[ii] * covXf[jj] * temp_kk) * fskew + covXf[ii] * covXf[jj] * covXf[kk] * temp_fskew -
1156 3.0 * covXf[ii] * covXf[jj] * covXf[kk] * fskew * temp_fvar / fvar) / fvar3;
1157 }
1158 }
1159 } // loop kk
1160 } // loop jj
1161 } // loop ii
1162
1163 rCM3vec[0] /= N;
1164
1165 UNPROTECT(1);
1166 return CM3vec;
1167 }
1168
CM3_CC(SEXP XXc,SEXP XXc2,SEXP mmargvars,SEXP mmargskews,SEXP mmargkurts,SEXP mmarg5s,SEXP mmarg6s,SEXP mm11,SEXP mm21,SEXP mm31,SEXP mm32,SEXP mm41,SEXP mm61,SEXP rr2,SEXP rr4,SEXP rr5,SEXP NN,SEXP PP)1169 SEXP CM3_CC(SEXP XXc, SEXP XXc2, SEXP mmargvars, SEXP mmargskews,
1170 SEXP mmargkurts, SEXP mmarg5s, SEXP mmarg6s,
1171 SEXP mm11, SEXP mm21, SEXP mm31, SEXP mm32,
1172 SEXP mm41, SEXP mm61, SEXP rr2, SEXP rr4, SEXP rr5,
1173 SEXP NN, SEXP PP){
1174 /*
1175 arguments
1176 XXc : numeric vector with centered observations
1177 XXc2 : numeric vector with centered and squared observations
1178 mmargvars : numeric vector with marginal variances
1179 mmargskews : numeric vector with marginal skewness values
1180 mmargkurts : numeric vector with marginal kurtosis values
1181 mmarg5s : numeric vector with marginal 5th order moments
1182 mmarg6s : numeric vector with marginal 6th order moments
1183 mm11 : numeric vector of t(Xc) %*% Xc / NN
1184 mm21 : numeric vector of t(Xc^2) %*% Xc / NN
1185 mm31 : numeric vector of t(Xc^3) %*% Xc / NN
1186 mm32 : numeric vector of t(Xc^3) %*% Xc^2 / NN
1187 mm41 : numeric vector of t(Xc^4) %*% Xc / NN
1188 mm61 : numeric vector of t(Xc^6) %*% Xc / NN
1189 rr2 : generalized correlation coefficient of order 2
1190 rr4 : generalized correlation coefficient of order 4
1191 rr5 : generalized correlation coefficient of order 5
1192 NN : integer, number of observations
1193 PP : integer, number of assets
1194
1195 Written by Dries Cornilly
1196 */
1197
1198 // // declare pointers for the vector
1199 double *Xc, *Xc2, *m11, *m21, *m31, *m32, *m41, *m61, *margvars, *margskews, *margkurts, *marg5s, *marg6s;
1200
1201 // use REAL() to access the C array inside the numeric vector passed in from R
1202 Xc = REAL(XXc);
1203 Xc2 = REAL(XXc2);
1204 m11 = REAL(mm11);
1205 m21 = REAL(mm21);
1206 m31 = REAL(mm31);
1207 m32 = REAL(mm32);
1208 m41 = REAL(mm41);
1209 m61 = REAL(mm61);
1210 margvars = REAL(mmargvars);
1211 margskews = REAL(mmargskews);
1212 margkurts = REAL(mmargkurts);
1213 marg5s = REAL(mmarg5s);
1214 marg6s = REAL(mmarg6s);
1215 double r2 = asReal(rr2);
1216 double r4 = asReal(rr4);
1217 double r5 = asReal(rr5);
1218 double N = asReal(NN);
1219 int n = asInteger(NN);
1220 int P = asInteger(PP);
1221
1222 // allocate and compute the sample variance and covariances
1223 SEXP CM3vec = PROTECT(allocVector(REALSXP, 1));
1224 double *rCM3vec = REAL(CM3vec);
1225 rCM3vec[0] = 0.0;
1226
1227 // loop over the unique indices i <= j <= k
1228 for (int ii = 0; ii < P; ii++) {
1229 int iiP = ii * P;
1230 int iiN = ii * n;
1231 for (int jj = ii; jj < P; jj++) {
1232 int jjP = jj * P;
1233 int jjN = jj * n;
1234 for (int kk = jj; kk < P; kk++) {
1235 int kkP = kk * P;
1236 int kkN = kk * n;
1237 if (ii == jj) {
1238 if (jj == kk) {
1239 // element phi_iii
1240 rCM3vec[0] += (marg6s[ii] - margskews[ii] * margskews[ii] - 6.0 * margkurts[ii] * margvars[ii] +
1241 9.0 * margvars[ii] * margvars[ii] * margvars[ii]);
1242 } else {
1243 // element phi_iik
1244 double temp_ii4 = m61[kkP + ii] - m21[kkP + ii] * margkurts[ii] - 4.0 * m31[kkP + ii] * margskews[ii] -
1245 2.0 * m11[kkP + ii] * marg5s[ii] - margvars[ii] * m41[kkP + ii] +
1246 12.0 * margvars[ii] * m11[kkP + ii] * margskews[ii];
1247 double temp_kk2 = m32[iiP + kk] - margvars[kk] * m21[kkP + ii] - margskews[kk] * margvars[ii] -
1248 2.0 * m21[iiP + kk] * m11[kkP + ii];
1249
1250 rCM3vec[0] += 3.0 * r2 * (sqrt(margvars[kk] / margkurts[ii]) * temp_ii4 +
1251 sqrt(margkurts[ii] / margvars[kk]) * temp_kk2) / 2.0;
1252 }
1253 } else {
1254 if (jj == kk) {
1255 // element phi_ijj
1256 double temp_ii2 = m32[jjP + ii] - margvars[ii] * m21[iiP + jj] - margskews[ii] * margvars[jj] -
1257 2.0 * m21[jjP + ii] * m11[jjP + ii];
1258 double temp_jj4 = m61[iiP + jj] - m21[iiP + jj] * margkurts[jj] - 4.0 * m31[iiP + jj] * margskews[jj] -
1259 2.0 * m11[jjP + ii] * marg5s[jj] - margvars[jj] * m41[iiP + jj] +
1260 12.0 * margvars[jj] * m11[jjP + ii] * margskews[jj];
1261
1262 rCM3vec[0] += 3.0 * r2 * (sqrt(margvars[ii] / margkurts[jj]) * temp_jj4 +
1263 sqrt(margkurts[jj] / margvars[ii]) * temp_ii2) / 2.0;
1264 } else {
1265 // element phi_ijk
1266 double S511 = 0.0;
1267 double S151 = 0.0;
1268 double S115 = 0.0;
1269 double S211 = 0.0;
1270 double S121 = 0.0;
1271 double S112 = 0.0;
1272 double S311 = 0.0;
1273 double S131 = 0.0;
1274 double S113 = 0.0;
1275 double m111 = 0.0;
1276 for (int tt = 0; tt < n; tt++) {
1277 S511 += Xc2[iiN + tt] * Xc2[iiN + tt] * Xc[iiN + tt] * Xc[jjN + tt] * Xc[kkN + tt];
1278 S151 += Xc[iiN + tt] * Xc2[jjN + tt] * Xc2[jjN + tt] * Xc[jjN + tt] * Xc[kkN + tt];
1279 S115 += Xc[iiN + tt] * Xc[jjN + tt] * Xc2[kkN + tt] * Xc2[kkN + tt] * Xc[kkN + tt];
1280 S211 += Xc2[iiN + tt] * Xc[jjN + tt] * Xc[kkN + tt];
1281 S121 += Xc[iiN + tt] * Xc2[jjN + tt] * Xc[kkN + tt];
1282 S112 += Xc[iiN + tt] * Xc[jjN + tt] * Xc2[kkN + tt];
1283 S311 += Xc2[iiN + tt] * Xc[iiN + tt] * Xc[jjN + tt] * Xc[kkN + tt];
1284 S131 += Xc[iiN + tt] * Xc2[jjN + tt] * Xc[jjN + tt] * Xc[kkN + tt];
1285 S113 += Xc[iiN + tt] * Xc[jjN + tt] * Xc2[kkN + tt] * Xc[kkN + tt];
1286 m111 += Xc[iiN + tt] * Xc[jjN + tt] * Xc[kkN + tt];
1287 }
1288 m111 /= N;
1289
1290 double temp_ii4 = S511 / N - m111 * margkurts[ii] - 4.0 * margskews[ii] * S211 / N -
1291 marg5s[ii] * m11[kkP + jj] - m41[jjP + ii] * m11[kkP + ii] - m41[kkP + ii] * m11[jjP + ii] +
1292 8.0 * margskews[ii] * m11[jjP + ii] * m11[kkP + ii] + 4.0 * margskews[ii] * margvars[ii] * m11[kkP + jj];
1293 double temp_jj4 = S151 / N - m111 * margkurts[jj] - 4.0 * margskews[jj] * S121 / N -
1294 marg5s[jj] * m11[kkP + ii] - m41[kkP + jj] * m11[iiP + jj] - m41[iiP + jj] * m11[kkP + jj] +
1295 8.0 * margskews[jj] * m11[kkP + jj] * m11[iiP + jj] + 4.0 * margskews[jj] * margvars[jj] * m11[kkP + ii];
1296 double temp_kk4 = S115 / N - m111 * margkurts[kk] - 4.0 * margskews[kk] * S112 / N -
1297 marg5s[kk] * m11[jjP + ii] - m41[iiP + kk] * m11[jjP + kk] - m41[jjP + kk] * m11[iiP + kk] +
1298 8.0 * margskews[kk] * m11[iiP + kk] * m11[jjP + kk] + 4.0 * margskews[kk] * margvars[kk] * m11[jjP + ii];
1299
1300 double temp_ii2 = S311 / N - margvars[ii] * m111 - margskews[ii] * m11[kkP + jj] -
1301 m21[jjP + ii] * m11[kkP + ii] - m21[kkP + ii] * m11[jjP + ii];
1302 double temp_jj2 = S131 / N - margvars[jj] * m111 - margskews[jj] * m11[kkP + ii] -
1303 m21[iiP + jj] * m11[kkP + jj] - m21[kkP + jj] * m11[iiP + jj];
1304 double temp_kk2 = S113 / N - margvars[kk] * m111 - margskews[kk] * m11[jjP + ii] -
1305 m21[iiP + kk] * m11[jjP + kk] - m21[jjP + kk] * m11[iiP + kk];
1306
1307 double cov_ii = sqrt(margvars[ii] * sqrt(margkurts[jj] / (margkurts[kk] * margkurts[kk] * margkurts[kk]))) * temp_kk4 / 2.0 +
1308 sqrt(margvars[ii] * sqrt(margkurts[kk] / (margkurts[jj] * margkurts[jj] * margkurts[jj]))) * temp_jj4 / 2.0 +
1309 sqrt(sqrt(margkurts[jj] * margkurts[kk]) / margvars[ii]) * temp_ii2;
1310 double cov_jj = sqrt(margvars[jj] * sqrt(margkurts[ii] / (margkurts[kk] * margkurts[kk] * margkurts[kk]))) * temp_kk4 / 2.0 +
1311 sqrt(margvars[jj] * sqrt(margkurts[kk] / (margkurts[ii] * margkurts[ii] * margkurts[ii]))) * temp_ii4 / 2.0 +
1312 sqrt(sqrt(margkurts[ii] * margkurts[kk]) / margvars[jj]) * temp_jj2;
1313 double cov_kk = sqrt(margvars[kk] * sqrt(margkurts[jj] / (margkurts[ii] * margkurts[ii] * margkurts[ii]))) * temp_ii4 / 2.0 +
1314 sqrt(margvars[kk] * sqrt(margkurts[ii] / (margkurts[jj] * margkurts[jj] * margkurts[jj]))) * temp_jj4 / 2.0 +
1315 sqrt(sqrt(margkurts[ii] * margkurts[jj]) / margvars[kk]) * temp_kk2;
1316
1317 rCM3vec[0] += r4 * sqrt(r5) * (cov_ii + cov_jj + cov_kk);
1318 }
1319 }
1320 } // loop kk
1321 } // loop jj
1322 } // loop ii
1323
1324 rCM3vec[0] /= N;
1325
1326 UNPROTECT(1);
1327 return CM3vec;
1328 }
1329
1330
1331 // // //
1332 // // Shrinkage and structured estimators for cokurtosis
1333 // // //
1334
M4_T12(SEXP mmargk_iiii,SEXP mmargk_iikk,SEXP PP)1335 SEXP M4_T12(SEXP mmargk_iiii, SEXP mmargk_iikk, SEXP PP){
1336 /*
1337 arguments
1338 mmargk_iiii : vector of length PP with marginal kurtosis values
1339 mmargk_iikk : vector of length PP with k_iikk values
1340 PP : integer, number of assets
1341
1342 Written by Dries Cornilly
1343 */
1344
1345 // // declare pointers for the vector
1346 double *margk_iiii, *margk_iikk;
1347
1348 // use REAL() to access the C array inside the numeric vector passed in from R
1349 margk_iiii = REAL(mmargk_iiii);
1350 margk_iikk = REAL(mmargk_iikk);
1351 int P = asInteger(PP);
1352
1353 // allocate and compute the coskewness matrix
1354 SEXP M4 = PROTECT(allocVector(REALSXP, P * (P + 1) * (P + 2) * (P + 3) / 24));
1355 double *rM4 = REAL(M4);
1356
1357 int iter = 0;
1358 // loop over the unique indices i <= j <= k <= l
1359 for (int ii = 0; ii < P; ii++) {
1360 for (int jj = ii; jj < P; jj++) {
1361 for (int kk = jj; kk < P; kk++) {
1362 for (int ll = kk; ll < P; ll++) {
1363 // compute coskewness element
1364 double elem = 0.0;
1365 if (((ii == jj) && (jj == kk)) && (kk == ll)) elem = margk_iiii[ii];
1366 if (((ii == jj) && (kk == ll)) && (jj != kk)) elem = margk_iikk[ii] * margk_iikk[kk];
1367
1368 // fill in the element and update the iterator
1369 rM4[iter] = elem;
1370 iter++;
1371 } // loop ll
1372 } // loop kk
1373 } // loop jj
1374 } // loop ii
1375
1376 UNPROTECT(1);
1377 return M4;
1378 }
1379
M4_CCoefficients(SEXP mmargvars,SEXP mmargkurts,SEXP mmarg6s,SEXP mm22,SEXP mm31,SEXP XXc,SEXP NN,SEXP PP)1380 SEXP M4_CCoefficients(SEXP mmargvars, SEXP mmargkurts, SEXP mmarg6s,
1381 SEXP mm22, SEXP mm31, SEXP XXc, SEXP NN, SEXP PP){
1382 /*
1383 arguments
1384 mmargvars : vector of length PP with marginal variances
1385 mmargkurts : vector of length PP with marginal kurtosis values
1386 mmarg6s : vector of length PP with marginal 6th order moments
1387 mm22 : t(Xc^2) %*% Xc^2
1388 mm31 : t(Xc^3) %*% Xc
1389 XXc : numeric vector with centered observations
1390 NN : integer, number of observations
1391 PP : integer, number of assets
1392
1393 Written by Dries Cornilly
1394 */
1395
1396 // // declare pointers for the vector
1397 double *margvars, *margkurts, *marg6s, *m22, *m31, *Xc;
1398
1399 // use REAL() to access the C array inside the numeric vector passed in from R
1400 margvars = REAL(mmargvars);
1401 margkurts = REAL(mmargkurts);
1402 marg6s = REAL(mmarg6s);
1403 m22 = REAL(mm22);
1404 m31 = REAL(mm31);
1405 Xc = REAL(XXc);
1406 double N = asReal(NN);
1407 int n = asInteger(NN);
1408 int P = asInteger(PP);
1409 double p = asReal(PP);
1410
1411 // allocate and compute the generalized correlation coefficients
1412 SEXP CCoef = PROTECT(allocVector(REALSXP, 4));
1413 double *rCCoef = REAL(CCoef);
1414
1415 // compute the generalized correlation coefficients
1416 rCCoef[0] = 0.0;
1417 rCCoef[1] = 0.0;
1418 for (int ii = 0; ii < P; ii++) {
1419 for (int jj = ii + 1; jj < P; jj++) {
1420 int jjP = jj * P;
1421 rCCoef[0] += m31[jjP + ii] / sqrt(marg6s[ii] * margvars[jj]);
1422 rCCoef[1] += m22[jjP + ii] / sqrt(margkurts[ii] * margkurts[jj]);
1423 }
1424 }
1425 rCCoef[0] *= 2.0 / (p * (p - 1.0));
1426 rCCoef[1] *= 2.0 / (p * (p - 1.0));
1427
1428 rCCoef[2] = 0.0;
1429 rCCoef[3] = 0.0;
1430 for (int ii = 0; ii < P; ii++) {
1431 int iiN = ii * n;
1432 for (int jj = ii + 1; jj < P; jj++) {
1433 int jjN = jj * n;
1434 for (int kk = jj + 1; kk < P; kk++) {
1435 int kkN = kk * n;
1436 double m211 = 0.0;
1437 for (int tt = 0; tt < n; tt++) {
1438 m211 += Xc[iiN + tt] * Xc[iiN + tt] * Xc[jjN + tt] * Xc[kkN + tt];
1439 }
1440 m211 /= N;
1441 rCCoef[2] += m211 / sqrt(margkurts[ii] * rCCoef[1] * sqrt(margkurts[jj] * margkurts[kk]));
1442
1443 for (int ll = kk + 1; ll < P; ll++) {
1444 int llN = ll * n;
1445 double m1111 = 0.0;
1446 for (int tt = 0; tt < n; tt++) {
1447 m1111 += Xc[iiN + tt] * Xc[jjN + tt] * Xc[kkN + tt] * Xc[llN + tt];
1448 }
1449 m1111 /= N;
1450 rCCoef[3] += m1111 / (rCCoef[1] * sqrt(sqrt(margkurts[ii] * margkurts[jj] * margkurts[kk] * margkurts[ll])));
1451 }
1452 }
1453 }
1454 }
1455 rCCoef[2] *= 6.0 / (p * (p - 1.0) * (p - 2.0));
1456 rCCoef[3] *= 24.0 / (p * (p - 1.0) * (p - 2.0) * (p - 3.0));
1457
1458 UNPROTECT(1);
1459 return CCoef;
1460 }
1461
M4_CC(SEXP mmargvars,SEXP mmargkurts,SEXP mmarg6s,SEXP rr3,SEXP rr5,SEXP rr6,SEXP rr7,SEXP PP)1462 SEXP M4_CC(SEXP mmargvars, SEXP mmargkurts, SEXP mmarg6s,
1463 SEXP rr3, SEXP rr5, SEXP rr6, SEXP rr7, SEXP PP){
1464 /*
1465 arguments
1466 mmargvars : vector of length PP with marginal variances
1467 mmargkurts : vector of length PP with marginal kurtosis values
1468 mmarg6s : vector of length PP with marginal 6th order central moments
1469 rr3 : generalized correlation of order 3
1470 rr5 : generalized correlation of order 5
1471 rr6 : generalized correlation of order 6
1472 rr7 : generalized correlation of order 7
1473 PP : integer, number of assets
1474
1475 Written by Dries Cornilly
1476 */
1477
1478 // // declare pointers for the vector
1479 double *margvars, *margkurts, *marg6s;
1480
1481 // use REAL() to access the C array inside the numeric vector passed in from R
1482 margvars = REAL(mmargvars);
1483 margkurts = REAL(mmargkurts);
1484 marg6s = REAL(mmarg6s);
1485 double r3 = asReal(rr3);
1486 double r5 = asReal(rr5);
1487 double r6 = asReal(rr6);
1488 double r7 = asReal(rr7);
1489 int P = asInteger(PP);
1490
1491 // allocate and compute the cokurtosis matrix
1492 SEXP M4 = PROTECT(allocVector(REALSXP, P * (P + 1) * (P + 2) * (P + 3) / 24));
1493 double *rM4 = REAL(M4);
1494
1495 int iter = 0;
1496 // loop over the unique indices i <= j <= k <= l
1497 for (int ii = 0; ii < P; ii++) {
1498 for (int jj = ii; jj < P; jj++) {
1499 for (int kk = jj; kk < P; kk++) {
1500 for (int ll = kk; ll < P; ll++) {
1501 // compute and fill cokurtosis element
1502 if (ii == jj) {
1503 if (jj == kk) {
1504 if (kk == ll) {
1505 // psi_iiii
1506 rM4[iter] = margkurts[ii];
1507 } else {
1508 // psi_iiil
1509 rM4[iter] = r3 * sqrt(marg6s[ii] * margvars[ll]);
1510 }
1511 } else {
1512 if (kk == ll) {
1513 // psi_iikk
1514 rM4[iter] = r5 * sqrt(margkurts[ii] * margkurts[kk]);
1515 } else {
1516 // psi_iikl
1517 rM4[iter] = r6 * sqrt(margvars[ii] * r5 * sqrt(margkurts[kk] * margkurts[ll]));
1518 }
1519 }
1520 } else {
1521 if (jj == kk) {
1522 if (kk == ll) {
1523 // psi_ijjj
1524 rM4[iter] = r3 * sqrt(margvars[ii] * marg6s[jj]);
1525 } else {
1526 // psi_ijjl
1527 rM4[iter] = r6 * sqrt(margvars[jj] * r5 * sqrt(margkurts[ii] * margkurts[ll]));
1528 }
1529 } else {
1530 if (kk == ll) {
1531 // psi_ijkk
1532 rM4[iter] = r6 * sqrt(margvars[kk] * r5 * sqrt(margkurts[ii] * margkurts[jj]));
1533 } else {
1534 // psi_ijkl
1535 rM4[iter] = r7 * r5 * sqrt(sqrt(margkurts[ii] * margkurts[jj] * margkurts[kk] * margkurts[ll]));
1536 }
1537 }
1538 }
1539 iter++;
1540 } // loop ii
1541 } // loop kk
1542 } // loop jj
1543 } // loop ii
1544
1545 UNPROTECT(1);
1546 return M4;
1547 }
1548
M4_1f(SEXP mmargkurts,SEXP ffvar,SEXP ffkurt,SEXP eepsvars,SEXP bbeta,SEXP PP)1549 SEXP M4_1f(SEXP mmargkurts, SEXP ffvar, SEXP ffkurt, SEXP eepsvars, SEXP bbeta, SEXP PP){
1550 /*
1551 arguments
1552 mmargkurts : vector of length PP with marginal kurtosis values
1553 ffvar : double with the factor variance
1554 ffkurt : vector with the factor fourth order central moment
1555 eepsvars : numeric vector with the variances of the idiosyncratic term
1556 bbeta : regression coefficients (factor loadings)
1557 PP : integer, number of assets
1558
1559 Written by Dries Cornilly
1560 */
1561
1562 // // declare pointers for the vector
1563 double *margkurts, *epsvars, *beta;
1564
1565 // use REAL() to access the C array inside the numeric vector passed in from R
1566 margkurts = REAL(mmargkurts);
1567 epsvars = REAL(eepsvars);
1568 beta = REAL(bbeta);
1569 double fvar = asReal(ffvar);
1570 double fkurt = asReal(ffkurt);
1571 int P = asInteger(PP);
1572
1573 // allocate and compute the cokurtosis matrix
1574 SEXP M4 = PROTECT(allocVector(REALSXP, P * (P + 1) * (P + 2) * (P + 3) / 24));
1575 double *rM4 = REAL(M4);
1576
1577 int iter = 0;
1578 // loop over the unique indices i <= j <= k <= l
1579 for (int ii = 0; ii < P; ii++) {
1580 for (int jj = ii; jj < P; jj++) {
1581 for (int kk = jj; kk < P; kk++) {
1582 for (int ll = kk; ll < P; ll++) {
1583 // compute and fill cokurtosis element
1584 if (ii == jj) {
1585 if (jj == kk) {
1586 if (kk == ll) {
1587 // psi_iiii
1588 rM4[iter] = margkurts[ii];
1589 } else {
1590 // psi_iiil
1591 rM4[iter] = beta[ii] * beta[ii] * beta[ii] * beta[ll] * fkurt +
1592 3.0 * beta[ii] * beta[ll] * fvar * epsvars[ii];
1593 }
1594 } else {
1595 if (kk == ll) {
1596 // psi_iikk
1597 rM4[iter] = beta[ii] * beta[ii] * beta[kk] * beta[kk] * fkurt +
1598 fvar * (beta[ii] * beta[ii] * epsvars[kk] + beta[kk] * beta[kk] * epsvars[ii]) +
1599 epsvars[ii] * epsvars[kk];
1600 } else {
1601 // psi_iikl
1602 rM4[iter] = beta[ii] * beta[ii] * beta[kk] * beta[ll] * fkurt +
1603 beta[kk] * beta[ll] * fvar * epsvars[ii];
1604 }
1605 }
1606 } else {
1607 if (jj == kk) {
1608 if (kk == ll) {
1609 // psi_ijjj
1610 rM4[iter] = beta[ii] * beta[jj] * beta[jj] * beta[jj] * fkurt +
1611 3.0 * beta[ii] * beta[jj] * fvar * epsvars[jj];
1612 } else {
1613 // psi_ijjl
1614 rM4[iter] = beta[ii] * beta[jj] * beta[jj] * beta[ll] * fkurt +
1615 beta[ii] * beta[ll] * fvar * epsvars[jj];
1616 }
1617 } else {
1618 if (kk == ll) {
1619 // psi_ijkk
1620 rM4[iter] = beta[ii] * beta[jj] * beta[kk] * beta[kk] * fkurt +
1621 beta[ii] * beta[jj] * fvar * epsvars[kk];
1622 } else {
1623 // psi_ijkl
1624 rM4[iter] = beta[ii] * beta[jj] * beta[kk] * beta[ll] * fkurt;
1625 }
1626 }
1627 }
1628 iter++;
1629 } // loop ii
1630 } // loop kk
1631 } // loop jj
1632 } // loop ii
1633
1634 UNPROTECT(1);
1635 return M4;
1636 }
1637
M4_MFresid(SEXP SStransf,SEXP eepsvars,SEXP PP)1638 SEXP M4_MFresid(SEXP SStransf, SEXP eepsvars, SEXP PP){
1639 /*
1640 arguments
1641 SStransf : contains the vectorized matrix B %*% S_F %*% t(B)
1642 eepsvars : numeric vector with the variances of the idiosyncratic term
1643 PP : integer, number of assets
1644
1645 Written by Dries Cornilly
1646 */
1647
1648 // // declare pointers for the vector
1649 double *Stransf, *epsvars;
1650
1651 // use REAL() to access the C array inside the numeric vector passed in from R
1652 epsvars = REAL(eepsvars);
1653 Stransf = REAL(SStransf);
1654 int P = asInteger(PP);
1655
1656 // allocate and compute the cokurtosis matrix
1657 SEXP M4 = PROTECT(allocVector(REALSXP, P * (P + 1) * (P + 2) * (P + 3) / 24));
1658 double *rM4 = REAL(M4);
1659
1660 int iter = 0;
1661 // loop over the unique indices i <= j <= k <= l
1662 for (int ii = 0; ii < P; ii++) {
1663 int iiP = ii * P;
1664 for (int jj = ii; jj < P; jj++) {
1665 for (int kk = jj; kk < P; kk++) {
1666 int kkP = kk * P;
1667 for (int ll = kk; ll < P; ll++) {
1668 // compute and fill cokurtosis element
1669 if (ii == jj) {
1670 if (jj == kk) {
1671 if (kk == ll) {
1672 // psi_iiii
1673 rM4[iter] = 6.0 * Stransf[iiP + ii] * epsvars[ii];
1674 } else {
1675 // psi_iiil
1676 rM4[iter] = 3.0 * Stransf[iiP + ll] * epsvars[ii];
1677 }
1678 } else {
1679 if (kk == ll) {
1680 // psi_iikk
1681 rM4[iter] = Stransf[iiP + ii] * epsvars[kk] + Stransf[kkP + kk] * epsvars[ii];
1682 } else {
1683 // psi_iikl
1684 rM4[iter] = Stransf[kkP + ll] * epsvars[ii];
1685 }
1686 }
1687 } else {
1688 if (jj == kk) {
1689 if (kk == ll) {
1690 // psi_ijjj
1691 rM4[iter] = 3.0 * Stransf[iiP + jj] * epsvars[jj];
1692 } else {
1693 // psi_ijjl
1694 rM4[iter] = Stransf[iiP + ll] * epsvars[jj];
1695 }
1696 } else {
1697 if (kk == ll) {
1698 // psi_ijkk
1699 rM4[iter] = Stransf[iiP + jj] * epsvars[kk];
1700 } else {
1701 // psi_ijkl
1702 rM4[iter] = 0.0;
1703 }
1704 }
1705 }
1706 iter++;
1707 } // loop ii
1708 } // loop kk
1709 } // loop jj
1710 } // loop ii
1711
1712 UNPROTECT(1);
1713 return M4;
1714 }
1715
VM4(SEXP XXc,SEXP XXc2,SEXP mm11,SEXP mm21,SEXP mm22,SEXP mm31,SEXP mm32,SEXP mm41,SEXP mm42,SEXP NN,SEXP PP)1716 SEXP VM4(SEXP XXc, SEXP XXc2, SEXP mm11, SEXP mm21, SEXP mm22, SEXP mm31,
1717 SEXP mm32, SEXP mm41, SEXP mm42, SEXP NN, SEXP PP){
1718 /*
1719 arguments
1720 XXc : numeric vector with centered observations
1721 XXc2 : numeric vector with centered and squared observations
1722 mm11 : numeric vector of t(Xc) %*% Xc / NN
1723 mm21 : numeric vector of t(Xc^2) %*% Xc / NN
1724 mm22 : numeric vector of t(Xc^2) %*% Xc^2 / NN
1725 mm31 : numeric vector of t(Xc^3) %*% Xc / NN
1726 mm32 : numeric vector of t(Xc^3) %*% Xc^2 / NN
1727 mm41 : numeric vector of t(Xc^4) %*% Xc / NN
1728 mm42 : numeric vector of t(Xc^4) %*% Xc^2 / NN
1729 NN : integer, number of observations
1730 PP : integer, number of assets
1731
1732 Written by Dries Cornilly
1733 */
1734
1735 // // declare pointers for the vector
1736 double *Xc, *Xc2, *m11, *m21, *m22, *m31, *m32, *m41, *m42;
1737
1738 // use REAL() to access the C array inside the numeric vector passed in from R
1739 Xc = REAL(XXc);
1740 Xc2 = REAL(XXc2);
1741 m11 = REAL(mm11);
1742 m21 = REAL(mm21);
1743 m22 = REAL(mm22);
1744 m31 = REAL(mm31);
1745 m32 = REAL(mm32);
1746 m41 = REAL(mm41);
1747 m42 = REAL(mm42);
1748 double N = asReal(NN);
1749 int n = asInteger(NN);
1750 int P = asInteger(PP);
1751 double p = asReal(PP);
1752
1753 // allocate and compute the sample variance and covariances
1754 SEXP VM4vec = PROTECT(allocVector(REALSXP, 3));
1755 double *rVM4vec = REAL(VM4vec);
1756 for (int ii = 0; ii < 3; ii++) {
1757 rVM4vec[ii] = 0.0;
1758 }
1759
1760 // loop over the unique indices i <= j <= k <= l
1761 for (int ii = 0; ii < P; ii++) {
1762 int iiP = ii * P;
1763 int iiN = ii * n;
1764 for (int jj = ii; jj < P; jj++) {
1765 int jjP = jj * P;
1766 int jjN = jj * n;
1767 for (int kk = jj; kk < P; kk++) {
1768 int kkP = kk * P;
1769 int kkN = kk * n;
1770 for (int ll = kk; ll < P; ll++) {
1771 int llP = ll * P;
1772 int llN = ll * n;
1773 if (ii == jj) {
1774 if (jj == kk) {
1775 if (kk == ll) {
1776 // psi_iiii
1777 double S8 = 0.0;
1778 for (int tt = 0; tt < n; tt++) {
1779 S8 += Xc2[iiN + tt] * Xc2[iiN + tt] * Xc2[iiN + tt] * Xc2[iiN + tt];
1780 }
1781
1782 double temp = (S8 / N - m31[iiP + ii] * m31[iiP + ii] -
1783 8.0 * m32[iiP + ii] * m21[iiP + ii] +
1784 16.0 * m11[iiP + ii] * m21[iiP + ii] * m21[iiP + ii]) / N;
1785 rVM4vec[0] += temp;
1786 rVM4vec[1] += temp / p; // cov with T (equal marginals)
1787 rVM4vec[2] += temp; // cov with T(unequal marginals)
1788 } else {
1789 // psi_iiil
1790 double S62 = 0.0;
1791 for (int tt = 0; tt < n; tt++) {
1792 S62 += Xc2[iiN + tt] * Xc2[iiN + tt] * Xc2[iiN + tt] * Xc2[llN + tt];
1793 }
1794
1795 rVM4vec[0] += 4.0 * (S62 / N - m31[llP + ii] * m31[llP + ii] -
1796 6.0 * m41[llP + ii] * m21[llP + ii] - 2.0 * m32[llP + ii] * m21[iiP + ii] +
1797 6.0 * m11[llP + ii] * m21[llP + ii] * m21[iiP + ii] +
1798 6.0 * m11[iiP + ii] * m21[llP + ii] * m21[llP + ii] +
1799 m21[iiP + ii] * m21[iiP + ii] * m11[llP + ll] +
1800 3.0 * m21[llP + ii] * m21[llP + ii] * m11[iiP + ii]) / N;
1801 }
1802 } else {
1803 if (kk == ll) {
1804 // psi_iikk
1805 double S44 = 0.0;
1806 for (int tt = 0; tt < n; tt++) {
1807 S44 += Xc2[iiN + tt] * Xc2[iiN + tt] * Xc2[kkN + tt] * Xc2[kkN + tt];
1808 }
1809
1810 rVM4vec[0] += 6.0 * (S44 / N - m22[kkP + ii] * m22[kkP + ii] -
1811 4.0 * m32[kkP + ii] * m21[iiP + kk] - 4.0 * m32[iiP + kk] * m21[kkP + ii] +
1812 8.0 * m11[kkP + ii] * m21[iiP + kk] * m21[kkP + ii] +
1813 4.0 * m21[kkP + ii] * m21[kkP + ii] * m11[kkP + kk] +
1814 4.0 * m21[iiP + kk] * m21[iiP + kk] * m11[iiP + ii]) / N;
1815
1816 // cov with T (equal marginals)
1817 double temp = 0.0;
1818 for (int mm = 0; mm < P; mm++) {
1819 int mmP = mm * P;
1820 int mmN = mm * n;
1821 double S222 = 0.0;
1822 for (int tt = 0; tt < N; tt++) {
1823 S222 += Xc2[iiN + tt] * Xc2[kkN + tt] * Xc2[mmN + tt];
1824 }
1825 temp += 2.0 * m11[mmP + mm] * (S222 / N - m11[mmP + mm] * m22[kkP + ii] -
1826 2.0 * m21[iiP + mm] * m21[iiP + kk] - 2.0 * m21[kkP + mm] * m21[kkP + ii]) / N;
1827 }
1828 rVM4vec[1] += 6.0 * temp / p;
1829
1830 // cov with T (unequal marginals)
1831 double temp_ii = m42[kkP + ii] - m11[iiP + ii] * m22[kkP + ii] -
1832 2.0 * m21[iiP + ii] * m21[iiP + kk] - 2.0 * m21[kkP + ii] * m21[kkP + ii];
1833 double temp_kk = m42[iiP + kk] - m11[kkP + kk] * m22[iiP + kk] -
1834 2.0 * m21[kkP + kk] * m21[kkP + ii] - 2.0 * m21[iiP + kk] * m21[iiP + kk];
1835 rVM4vec[2] += 6.0 * (m11[kkP + kk] * temp_ii + m11[iiP + ii] * temp_kk) / N;
1836 } else {
1837 // psi_iikl
1838 double S422 = 0.0;
1839 double S311 = 0.0;
1840 double S221 = 0.0;
1841 double S212 = 0.0;
1842 double m211 = 0.0;
1843 double m111 = 0.0;
1844 for (int tt = 0; tt < n; tt++) {
1845 S422 += Xc2[iiN + tt] * Xc2[iiN + tt] * Xc2[kkN + tt] * Xc2[llN + tt];
1846 S311 += Xc2[iiN + tt] * Xc[iiN + tt] * Xc[kkN + tt] * Xc[llN + tt];
1847 S221 += Xc2[iiN + tt] * Xc2[kkN + tt] * Xc[llN + tt];
1848 S212 += Xc2[iiN + tt] * Xc[kkN + tt] * Xc2[llN + tt];
1849 m211 += Xc2[iiN + tt] * Xc[kkN + tt] * Xc[llN + tt];
1850 m111 += Xc[iiN + tt] * Xc[kkN + tt] * Xc[llN + tt];
1851 }
1852 m211 /= N;
1853 m111 /= N;
1854
1855 rVM4vec[0] += 12.0 * (S422 / N - m211 * m211 - 4.0 * S311 / N * m111 -
1856 2.0 * S221 / N * m21[llP + ii] - 2.0 * S212 / N * m21[kkP + ii] +
1857 4.0 * m11[llP + ii] * m21[kkP + ii] * m111 + 4.0 * m11[kkP + ii] * m111 * m21[llP + ii] +
1858 2.0 * m11[llP + kk] * m21[llP + ii] * m21[kkP + ii] +
1859 4.0 * m11[iiP + ii] * m111 * m111 + m21[kkP + ii] * m21[kkP + ii] * m11[llP + ll] +
1860 m21[llP + ii] * m21[llP + ii] * m11[kkP + kk]) / N;
1861 }
1862 }
1863 } else {
1864 if (jj == kk) {
1865 if (kk == ll) {
1866 // psi_ijjj
1867 double S62 = 0.0;
1868 for (int tt = 0; tt < n; tt++) {
1869 S62 += Xc2[jjN + tt] * Xc2[jjN + tt] * Xc2[jjN + tt] * Xc2[iiN + tt];
1870 }
1871
1872 rVM4vec[0] += 4.0 * (S62 / N - m31[iiP + jj] * m31[iiP + jj] -
1873 6.0 * m41[iiP + jj] * m21[iiP + jj] - 2.0 * m32[iiP + jj] * m21[jjP + jj] +
1874 6.0 * m11[iiP + jj] * m21[iiP + jj] * m21[jjP + jj] +
1875 6.0 * m11[jjP + jj] * m21[iiP + jj] * m21[iiP + jj] +
1876 m21[jjP + jj] * m21[jjP + jj] * m11[iiP + ii] +
1877 3.0 * m21[iiP + jj] * m21[iiP + jj] * m11[jjP + jj]) / N;
1878 } else {
1879 // psi_ijjl
1880 double S422 = 0.0;
1881 double S311 = 0.0;
1882 double S221 = 0.0;
1883 double S212 = 0.0;
1884 double m211 = 0.0;
1885 double m111 = 0.0;
1886 for (int tt = 0; tt < n; tt++) {
1887 S422 += Xc2[jjN + tt] * Xc2[jjN + tt] * Xc2[iiN + tt] * Xc2[llN + tt];
1888 S311 += Xc2[jjN + tt] * Xc[jjN + tt] * Xc[iiN + tt] * Xc[llN + tt];
1889 S221 += Xc2[jjN + tt] * Xc2[iiN + tt] * Xc[llN + tt];
1890 S212 += Xc2[jjN + tt] * Xc[iiN + tt] * Xc2[llN + tt];
1891 m211 += Xc2[jjN + tt] * Xc[iiN + tt] * Xc[llN + tt];
1892 m111 += Xc[jjN + tt] * Xc[iiN + tt] * Xc[llN + tt];
1893 }
1894 m211 /= N;
1895 m111 /= N;
1896
1897 rVM4vec[0] += 12.0 * (S422 / N - m211 * m211 - 4.0 * S311 / N * m111 -
1898 2.0 * S221 / N * m21[llP + jj] - 2.0 * S212 / N * m21[iiP + jj] +
1899 4.0 * m11[llP + jj] * m21[iiP + jj] * m111 + 4.0 * m11[iiP + jj] * m111 * m21[llP + jj] +
1900 2.0 * m11[llP + ii] * m21[llP + jj] * m21[iiP + jj] +
1901 4.0 * m11[jjP + jj] * m111 * m111 + m21[iiP + jj] * m21[iiP + jj] * m11[llP + ll] +
1902 m21[llP + jj] * m21[llP + jj] * m11[iiP + ii]) / N;
1903 }
1904 } else {
1905 if (kk == ll) {
1906 // psi_ijkk
1907 double S422 = 0.0;
1908 double S311 = 0.0;
1909 double S221 = 0.0;
1910 double S212 = 0.0;
1911 double m211 = 0.0;
1912 double m111 = 0.0;
1913 for (int tt = 0; tt < n; tt++) {
1914 S422 += Xc2[kkN + tt] * Xc2[kkN + tt] * Xc2[iiN + tt] * Xc2[jjN + tt];
1915 S311 += Xc2[kkN + tt] * Xc[kkN + tt] * Xc[iiN + tt] * Xc[jjN + tt];
1916 S221 += Xc2[kkN + tt] * Xc2[iiN + tt] * Xc[jjN + tt];
1917 S212 += Xc2[kkN + tt] * Xc[iiN + tt] * Xc2[jjN + tt];
1918 m211 += Xc2[kkN + tt] * Xc[iiN + tt] * Xc[jjN + tt];
1919 m111 += Xc[kkN + tt] * Xc[iiN + tt] * Xc[jjN + tt];
1920 }
1921 m211 /= N;
1922 m111 /= N;
1923
1924 rVM4vec[0] += 12.0 * (S422 / N - m211 * m211 - 4.0 * S311 / N * m111 -
1925 2.0 * S221 / N * m21[jjP + kk] - 2.0 * S212 / N * m21[iiP + kk] +
1926 4.0 * m11[jjP + kk] * m21[iiP + kk] * m111 + 4.0 * m11[iiP + kk] * m111 * m21[jjP + kk] +
1927 2.0 * m11[jjP + ii] * m21[jjP + kk] * m21[iiP + kk] +
1928 4.0 * m11[kkP + kk] * m111 * m111 + m21[iiP + kk] * m21[iiP + kk] * m11[jjP + jj] +
1929 m21[jjP + kk] * m21[jjP + kk] * m11[iiP + ii]) / N;
1930 } else {
1931 // psi_ijkl
1932 double S2222 = 0.0;
1933 double S2111 = 0.0;
1934 double S1211 = 0.0;
1935 double S1121 = 0.0;
1936 double S1112 = 0.0;
1937 double m1111 = 0.0;
1938 double m0111 = 0.0;
1939 double m1011 = 0.0;
1940 double m1101 = 0.0;
1941 double m1110 = 0.0;
1942 for (int tt = 0; tt < n; tt++) {
1943 S2222 += Xc2[iiN + tt] * Xc2[jjN + tt] * Xc2[kkN + tt] * Xc2[llN + tt];
1944 S2111 += Xc2[iiN + tt] * Xc[jjN + tt] * Xc[kkN + tt] * Xc[llN + tt];
1945 S1211 += Xc[iiN + tt] * Xc2[jjN + tt] * Xc[kkN + tt] * Xc[llN + tt];
1946 S1121 += Xc[iiN + tt] * Xc[jjN + tt] * Xc2[kkN + tt] * Xc[llN + tt];
1947 S1112 += Xc[iiN + tt] * Xc[jjN + tt] * Xc[kkN + tt] * Xc2[llN + tt];
1948 m1111 += Xc[iiN + tt] * Xc[jjN + tt] * Xc[kkN + tt] * Xc[llN + tt];
1949 m0111 += Xc[jjN + tt] * Xc[kkN + tt] * Xc[llN + tt];
1950 m1011 += Xc[iiN + tt] * Xc[kkN + tt] * Xc[llN + tt];
1951 m1101 += Xc[iiN + tt] * Xc[jjN + tt] * Xc[llN + tt];
1952 m1110 += Xc[iiN + tt] * Xc[jjN + tt] * Xc[kkN + tt];
1953 }
1954 m1111 /= N;
1955 m0111 /= N;
1956 m1011 /= N;
1957 m1101 /= N;
1958 m1110 /= N;
1959
1960 rVM4vec[0] += 24.0 * (S2222 / N - m1111 * m1111 - 2.0 * S2111 / N * m0111 -
1961 2.0 * S1211 / N * m1011 - 2.0 * S1121 / N * m1101 - 2.0 * S1112 / N * m1110 +
1962 2.0 * m11[llP + ii] * m1110 * m0111 + 2.0 * m11[llP + jj] * m1011 * m1110 +
1963 2.0 * m11[llP + kk] * m1101 * m1110 + 2.0 * m11[kkP + ii] * m0111 * m1101 +
1964 2.0 * m11[kkP + jj] * m1011 * m1101 + 2.0 * m11[jjP + ii] * m0111 * m1011 +
1965 m1110 * m1110 * m11[llP + ll] + m1101 * m1101 * m11[kkP + kk] +
1966 m1011 * m1011 * m11[jjP + jj] + m0111 * m0111 * m11[iiP + ii]) / N;
1967 }
1968 }
1969 }
1970 } // loop ll
1971 } // loop kk
1972 } // loop jj
1973 } // loop ii
1974
1975 // cov with T (equal marginals)
1976 for (int ii = 0; ii < P; ii++) {
1977 int iiP = ii * P;
1978 int iiN = ii * n;
1979 for (int jj = ii + 1; jj < P; jj++) {
1980 int jjP = jj * P;
1981 int jjN = jj * n;
1982 double S44 = 0.0;
1983 for (int tt = 0; tt < n; tt++) {
1984 S44 += Xc2[iiN + tt] * Xc2[iiN + tt] * Xc2[jjN + tt] * Xc2[jjN + tt];
1985 }
1986 rVM4vec[1] += 2 * (S44 / N - m22[iiP + ii] * m22[jjP + jj] -
1987 4.0 * m41[jjP + ii] * m21[jjP + jj] - 4.0 * m41[iiP + jj] * m21[iiP + ii] +
1988 16.0 * m21[iiP + ii] * m21[jjP + jj] * m11[jjP + ii]) / (N * p);
1989 } // loop over jj
1990 } // loop over ii
1991
1992 UNPROTECT(1);
1993 return VM4vec;
1994 }
1995
CM4_CC(SEXP XXc,SEXP XXc2,SEXP mm11,SEXP mm21,SEXP mm22,SEXP mm31,SEXP mm32,SEXP mm33,SEXP mm41,SEXP rr3,SEXP rr5,SEXP rr6,SEXP rr7,SEXP mmarg6s,SEXP mmarg7s,SEXP NN,SEXP PP)1996 SEXP CM4_CC(SEXP XXc, SEXP XXc2, SEXP mm11, SEXP mm21, SEXP mm22,
1997 SEXP mm31, SEXP mm32, SEXP mm33, SEXP mm41, SEXP rr3,
1998 SEXP rr5, SEXP rr6, SEXP rr7, SEXP mmarg6s, SEXP mmarg7s,
1999 SEXP NN, SEXP PP){
2000 /*
2001 arguments
2002 XXc : numeric vector with centered observations
2003 XXc2 : numeric vector with centered and squared observations
2004 mm11 : numeric vector of t(Xc) %*% Xc / NN
2005 mm21 : numeric vector of t(Xc^2) %*% Xc / NN
2006 mm22 : numeric vector of t(Xc^2) %*% Xc^2 / NN
2007 mm31 : numeric vector of t(Xc^3) %*% Xc / NN
2008 mm32 : numeric vector of t(Xc^3) %*% Xc^2 / NN
2009 mm33 : numeric vector of t(Xc^3) %*% Xc^3 / NN
2010 mm41 : numeric vector of t(Xc^4) %*% Xc / NN
2011 rr3 : double with generalized correlation coefficient of order 3
2012 rr5 : double with generalized correlation coefficient of order 5
2013 rr6 : double with generalized correlation coefficient of order 6
2014 rr7 : double with generalized correlation coefficient of order 7
2015 mmarg6s : numeric vector with marginal 6th order central moments
2016 mmarg7s : numeric vector with marginal 7th order central moments
2017 NN : integer, number of observations
2018 PP : integer, number of assets
2019
2020 Written by Dries Cornilly
2021 */
2022
2023 // // declare pointers for the vector
2024 double *Xc, *Xc2, *m11, *m21, *m22, *m31, *m32, *m33, *m41, *marg6s, *marg7s;
2025
2026 // use REAL() to access the C array inside the numeric vector passed in from R
2027 Xc = REAL(XXc);
2028 Xc2 = REAL(XXc2);
2029 m11 = REAL(mm11);
2030 m21 = REAL(mm21);
2031 m22 = REAL(mm22);
2032 m31 = REAL(mm31);
2033 m32 = REAL(mm32);
2034 m33 = REAL(mm33);
2035 m41 = REAL(mm41);
2036 marg6s = REAL(mmarg6s);
2037 marg7s = REAL(mmarg7s);
2038 double N = asReal(NN);
2039 int n = asInteger(NN);
2040 int P = asInteger(PP);
2041 double r3 = asReal(rr3);
2042 double r5 = asReal(rr5);
2043 double r6 = asReal(rr6);
2044 double r7 = asReal(rr7);
2045
2046 // allocate and compute the covariance
2047 SEXP CM4 = PROTECT(allocVector(REALSXP, 1));
2048 double *rCM4 = REAL(CM4);
2049 rCM4[0] = 0.0;
2050
2051 // loop over the unique indices i <= j <= k <= l
2052 for (int ii = 0; ii < P; ii++) {
2053 int iiP = ii * P;
2054 int iiN = ii * n;
2055 for (int jj = ii; jj < P; jj++) {
2056 int jjP = jj * P;
2057 int jjN = jj * n;
2058 for (int kk = jj; kk < P; kk++) {
2059 int kkP = kk * P;
2060 int kkN = kk * n;
2061 for (int ll = kk; ll < P; ll++) {
2062 int llP = ll * P;
2063 int llN = ll * n;
2064 if (ii == jj) {
2065 if (jj == kk) {
2066 if (kk == ll) {
2067 // psi_iiii
2068 double S8 = 0.0;
2069 for (int tt = 0; tt < n; tt++) {
2070 S8 += Xc2[iiN + tt] * Xc2[iiN + tt] * Xc2[iiN + tt] * Xc2[iiN + tt];
2071 }
2072
2073 rCM4[0] += (S8 / N - m31[iiP + ii] * m31[iiP + ii] -
2074 8.0 * m41[iiP + ii] * m21[iiP + ii] +
2075 16.0 * m11[iiP + ii] * m21[iiP + ii] * m21[iiP + ii]) / N;
2076 } else {
2077 // psi_iiil
2078 double S91 = 0.0;
2079 double S62 = 0.0;
2080 for (int tt = 0; tt < n; tt++) {
2081 S91 += Xc2[iiN + tt] * Xc2[iiN + tt] * Xc2[iiN + tt] *
2082 Xc2[iiN + tt] * Xc[iiN + tt] * Xc[llN + tt];
2083 S62 += Xc2[iiN + tt] * Xc2[iiN + tt] * Xc2[iiN + tt] * Xc2[llN + tt];
2084 }
2085
2086 double temp_ii = S91 / N - m31[llP + ii] * marg6s[ii] -
2087 3.0 * m21[llP + ii] * marg7s[ii] - m21[iiP + ii] * S62 / N -
2088 6.0 * m41[iiP + ii] * m41[llP + ii] +
2089 18.0 * m41[iiP + ii] * m21[llP + ii] * m11[iiP + ii] +
2090 6.0 * m41[iiP + ii] * m21[iiP + ii] * m11[llP + ii];
2091 double temp_ll = m33[llP + ii] - m31[llP + ii] * m11[llP + ll] -
2092 3.0 * m21[iiP + ll] * m21[llP + ii] - m21[iiP + ii] * m21[llP + ll];
2093
2094 rCM4[0] += 2.0 * r3 * (sqrt(m11[llP + ll] / marg6s[ii]) * temp_ii +
2095 sqrt(marg6s[ii] / m11[llP + ll]) * temp_ll) / N;
2096 }
2097 } else {
2098 if (kk == ll) {
2099 // psi_iikk
2100 double S62 = 0.0;
2101 double S26 = 0.0;
2102 for (int tt = 0; tt < n; tt++) {
2103 S62 += Xc2[iiN + tt] * Xc2[iiN + tt] * Xc2[iiN + tt] * Xc2[kkN + tt];
2104 S26 += Xc2[iiN + tt] * Xc2[kkN + tt] * Xc2[kkN + tt] * Xc2[kkN + tt];
2105 }
2106
2107 double temp_ii = S62 / N - m22[kkP + ii] * m22[iiP + ii] -
2108 4.0 * m32[kkP + ii] * m21[iiP + ii] - 2.0 * m21[iiP + kk] * m32[iiP + ii] -
2109 2.0 * m21[kkP + ii] * m41[kkP + ii] + 8.0 * m21[kkP + ii] * m21[iiP + ii] * m11[kkP + ii] +
2110 8.0 * m21[iiP + kk] * m21[iiP + ii] * m11[iiP + ii];
2111 double temp_kk = S26 / N - m22[iiP + kk] * m22[kkP + kk] -
2112 4.0 * m32[iiP + kk] * m21[kkP + kk] - 2.0 * m21[kkP + ii] * m32[kkP + kk] -
2113 2.0 * m21[iiP + kk] * m41[iiP + kk] + 8.0 * m21[iiP + kk] * m21[kkP + kk] * m11[iiP + kk] +
2114 8.0 * m21[kkP + ii] * m21[kkP + kk] * m11[kkP + kk];
2115
2116 rCM4[0] += 3.0 * r5 * (sqrt(m22[kkP + kk] / m22[iiP + ii]) * temp_ii +
2117 sqrt(m22[iiP + ii] / m22[kkP + kk]) * temp_kk) / N;
2118 } else {
2119 // psi_iikl
2120 double S611 = 0.0;
2121 double S251 = 0.0;
2122 double S215 = 0.0;
2123 double S311 = 0.0;
2124 double S221 = 0.0;
2125 double S212 = 0.0;
2126 double m211 = 0.0;
2127 double m111 = 0.0;
2128 for (int tt = 0; tt < n; tt++) {
2129 S611 += Xc2[iiN + tt] * Xc2[iiN + tt] * Xc2[iiN + tt] * Xc[kkN + tt] * Xc[llN + tt];
2130 S251 += Xc2[iiN + tt] * Xc2[kkN + tt] * Xc2[kkN + tt] * Xc[kkN + tt] * Xc[llN + tt];
2131 S215 += Xc2[iiN + tt] * Xc[kkN + tt] * Xc2[llN + tt] * Xc2[llN + tt] * Xc[llN + tt];
2132 S311 += Xc2[iiN + tt] * Xc[iiN + tt] * Xc[kkN + tt] * Xc[llN + tt];
2133 S221 += Xc2[iiN + tt] * Xc2[kkN + tt] * Xc[llN + tt];
2134 S212 += Xc2[iiN + tt] * Xc[kkN + tt] * Xc2[llN + tt];
2135 m211 += Xc2[iiN + tt] * Xc[kkN + tt] * Xc[llN + tt];
2136 m111 += Xc[iiN + tt] * Xc[kkN + tt] * Xc[llN + tt];
2137 }
2138 m211 /= N;
2139 m111 /= N;
2140
2141 double temp_ii = S611 / N - m211 * m22[iiP + ii] - 4.0 * S311 / N * m21[iiP + ii] -
2142 2.0 * m111 * m32[iiP + ii] - m21[llP + ii] * m41[kkP + ii] - m21[kkP + ii] * m41[llP + ii] +
2143 4.0 * m21[kkP + ii] * m21[iiP + ii] * m11[llP + ii] + 4.0 * m21[llP + ii] * m21[iiP + ii] * m11[kkP + ii] +
2144 8.0 * m111 * m21[iiP + ii] * m11[iiP + ii];
2145 double temp_kk = S251 / N - m211 * m22[kkP + kk] - 4.0 * S221 / N * m21[kkP + kk] -
2146 m21[llP + ii] * m32[kkP + kk] - 2.0 * m111 * m41[iiP + kk] - m21[kkP + ii] * m41[llP + kk] +
2147 4.0 * m21[kkP + ii] * m21[kkP + kk] * m11[llP + kk] + 8.0 * m111 * m21[kkP + kk] * m11[kkP + ii] +
2148 4.0 * m21[llP + ii] * m11[kkP + kk] * m21[kkP + kk];
2149 double temp_ll = S215 / N - m211 * m22[llP + ll] - 4.0 * S212 / N * m21[llP + ll] -
2150 m21[kkP + ii] * m32[llP + ll] - 2.0 * m111 * m41[iiP + ll] - m21[llP + ii] * m41[kkP + ll] +
2151 4.0 * m21[llP + ii] * m21[llP + ll] * m11[kkP + ll] + 8.0 * m111 * m21[llP + ll] * m11[llP + ii] +
2152 4.0 * m21[kkP + ii] * m11[llP + ll] * m21[llP + ll];
2153
2154 rCM4[0] += 3.0 * r6 * sqrt(r5) *
2155 (2.0 * sqrt(sqrt(m22[kkP + kk] * m22[llP + ll]) / m22[iiP + ii]) * temp_ii +
2156 sqrt(m22[iiP + ii] * sqrt(m22[llP + ll] / (m22[kkP + kk] * m22[kkP + kk] * m22[kkP + kk]))) * temp_kk +
2157 sqrt(m22[iiP + ii] * sqrt(m22[kkP + kk] / (m22[llP + ll] * m22[llP + ll] * m22[llP + ll]))) * temp_ll) / N;
2158 }
2159 }
2160 } else {
2161 if (jj == kk) {
2162 if (kk == ll) {
2163 // psi_ijjj
2164 double S91 = 0.0;
2165 double S62 = 0.0;
2166 for (int tt = 0; tt < n; tt++) {
2167 S91 += Xc2[jjN + tt] * Xc2[jjN + tt] * Xc2[jjN + tt] *
2168 Xc2[jjN + tt] * Xc[jjN + tt] * Xc[iiN + tt];
2169 S62 += Xc2[jjN + tt] * Xc2[jjN + tt] * Xc2[jjN + tt] * Xc2[iiN + tt];
2170 }
2171
2172 double temp_jj = S91 / N - m31[iiP + jj] * marg6s[jj] -
2173 3.0 * m21[iiP + jj] * marg7s[jj] - m21[jjP + jj] * S62 / N -
2174 6.0 * m41[jjP + jj] * m41[iiP + jj] +
2175 18.0 * m41[jjP + jj] * m21[iiP + jj] * m11[jjP + jj] +
2176 6.0 * m41[jjP + jj] * m21[jjP + jj] * m11[iiP + jj];
2177 double temp_ii = m33[iiP + jj] - m31[iiP + jj] * m11[iiP + ii] -
2178 3.0 * m21[jjP + ii] * m21[iiP + jj] - m21[jjP + jj] * m21[iiP + ii];
2179
2180 rCM4[0] += 2.0 * r3 * (sqrt(m11[iiP + ii] / marg6s[jj]) * temp_jj +
2181 sqrt(marg6s[jj] / m11[iiP + ii]) * temp_ii) / N;
2182 } else {
2183 // psi_ijjl
2184 double S161 = 0.0;
2185 double S521 = 0.0;
2186 double S125 = 0.0;
2187 double S131 = 0.0;
2188 double S221 = 0.0;
2189 double S122 = 0.0;
2190 double m121 = 0.0;
2191 double m111 = 0.0;
2192 for (int tt = 0; tt < n; tt++) {
2193 S161 += Xc[iiN + tt] * Xc2[jjN + tt] * Xc2[jjN + tt] * Xc2[jjN + tt] * Xc[llN + tt];
2194 S521 += Xc2[iiN + tt] * Xc2[iiN + tt] * Xc[iiN + tt] * Xc2[jjN + tt] * Xc[llN + tt];
2195 S125 += Xc[iiN + tt] * Xc2[jjN + tt] * Xc2[llN + tt] * Xc2[llN + tt] * Xc[llN + tt];
2196 S131 += Xc[iiN + tt] * Xc2[jjN + tt] * Xc[jjN + tt] * Xc[llN + tt];
2197 S221 += Xc2[iiN + tt] * Xc2[jjN + tt] * Xc[llN + tt];
2198 S122 += Xc[iiN + tt] * Xc2[jjN + tt] * Xc2[llN + tt];
2199 m121 += Xc[iiN + tt] * Xc2[jjN + tt] * Xc[llN + tt];
2200 m111 += Xc[iiN + tt] * Xc[jjN + tt] * Xc[llN + tt];
2201 }
2202 m121 /= N;
2203 m111 /= N;
2204
2205 double temp_jj = S161 / N - m121 * m22[jjP + jj] - 4.0 * S131 / N * m21[jjP + jj] -
2206 2.0 * m111 * m32[jjP + jj] - m21[llP + jj] * m41[iiP + jj] - m21[iiP + jj] * m41[llP + jj] +
2207 4.0 * m21[iiP + jj] * m21[jjP + jj] * m11[llP + jj] + 4.0 * m21[llP + jj] * m21[jjP + jj] * m11[iiP + jj] +
2208 8.0 * m111 * m21[jjP + jj] * m11[jjP + jj];
2209 double temp_ii = S521 / N - m121 * m22[iiP + ii] - 4.0 * S221 / N * m21[iiP + ii] -
2210 m21[llP + jj] * m32[iiP + ii] - 2.0 * m111 * m41[jjP + ii] - m21[iiP + jj] * m41[llP + ii] +
2211 4.0 * m21[iiP + jj] * m21[iiP + ii] * m11[llP + ii] + 8.0 * m111 * m21[iiP + ii] * m11[jjP + ii] +
2212 4.0 * m21[llP + jj] * m11[iiP + ii] * m21[iiP + ii];
2213 double temp_ll = S125 / N - m121 * m22[llP + ll] - 4.0 * S122 / N * m21[llP + ll] -
2214 m21[iiP + jj] * m32[llP + ll] - 2.0 * m111 * m41[jjP + ll] - m21[llP + jj] * m41[iiP + ll] +
2215 4.0 * m21[llP + jj] * m21[llP + ll] * m11[iiP + ll] + 8.0 * m111 * m21[llP + ll] * m11[jjP + ll] +
2216 4.0 * m21[iiP + jj] * m11[llP + ll] * m21[llP + ll];
2217
2218 rCM4[0] += 3.0 * r6 * sqrt(r5) *
2219 (2.0 * sqrt(sqrt(m22[iiP + ii] * m22[llP + ll]) / m22[jjP + jj]) * temp_jj +
2220 sqrt(m22[jjP + jj] * sqrt(m22[llP + ll] / (m22[iiP + ii] * m22[iiP + ii] * m22[iiP + ii]))) * temp_ii +
2221 sqrt(m22[jjP + jj] * sqrt(m22[iiP + ii] / (m22[llP + ll] * m22[llP + ll] * m22[llP + ll]))) * temp_ll) / N;
2222 }
2223 } else {
2224 if (kk == ll) {
2225 // psi_ijkk
2226 double S116 = 0.0;
2227 double S152 = 0.0;
2228 double S512 = 0.0;
2229 double S113 = 0.0;
2230 double S122 = 0.0;
2231 double S212 = 0.0;
2232 double m112 = 0.0;
2233 double m111 = 0.0;
2234 for (int tt = 0; tt < n; tt++) {
2235 S116 += Xc[iiN + tt] * Xc[jjN + tt] * Xc2[kkN + tt] * Xc2[kkN + tt] * Xc2[kkN + tt];
2236 S152 += Xc[iiN + tt] * Xc2[jjN + tt] * Xc2[jjN + tt] * Xc[jjN + tt] * Xc2[kkN + tt];
2237 S512 += Xc2[iiN + tt] * Xc2[iiN + tt] * Xc[iiN + tt] * Xc[jjN + tt] * Xc2[kkN + tt];
2238 S113 += Xc[iiN + tt] * Xc[jjN + tt] * Xc2[kkN + tt] * Xc[kkN + tt];
2239 S122 += Xc[iiN + tt] * Xc2[jjN + tt] * Xc2[kkN + tt];
2240 S212 += Xc2[iiN + tt] * Xc[jjN + tt] * Xc2[kkN + tt];
2241 m112 += Xc[iiN + tt] * Xc[jjN + tt] * Xc2[kkN + tt];
2242 m111 += Xc[iiN + tt] * Xc[jjN + tt] * Xc[kkN + tt];
2243 }
2244 m112 /= N;
2245 m111 /= N;
2246
2247 double temp_kk = S116 / N - m112 * m22[kkP + kk] - 4.0 * S113 / N * m21[kkP + kk] -
2248 2.0 * m111 * m32[kkP + kk] - m21[jjP + kk] * m41[iiP + kk] - m21[iiP + kk] * m41[jjP + kk] +
2249 4.0 * m21[iiP + kk] * m21[kkP + kk] * m11[jjP + kk] + 4.0 * m21[jjP + kk] * m21[kkP + kk] * m11[iiP + kk] +
2250 8.0 * m111 * m21[kkP + kk] * m11[kkP + kk];
2251 double temp_ii = S512 / N - m112 * m22[iiP + ii] - 4.0 * S212 / N * m21[iiP + ii] -
2252 m21[jjP + kk] * m32[iiP + ii] - 2.0 * m111 * m41[kkP + ii] - m21[iiP + kk] * m41[jjP + ii] +
2253 4.0 * m21[iiP + kk] * m21[iiP + ii] * m11[jjP + ii] + 8.0 * m111 * m21[iiP + ii] * m11[kkP + ii] +
2254 4.0 * m21[jjP + kk] * m11[iiP + ii] * m21[iiP + ii];
2255 double temp_jj = S152 / N - m112 * m22[jjP + jj] - 4.0 * S122 / N * m21[jjP + jj] -
2256 m21[iiP + kk] * m32[jjP + jj] - 2.0 * m111 * m41[kkP + jj] - m21[jjP + kk] * m41[iiP + jj] +
2257 4.0 * m21[jjP + kk] * m21[jjP + jj] * m11[iiP + jj] + 8.0 * m111 * m21[jjP + jj] * m11[kkP + jj] +
2258 4.0 * m21[iiP + kk] * m11[jjP + jj] * m21[jjP + jj];
2259
2260 rCM4[0] += 3.0 * r6 * sqrt(r5) *
2261 (2.0 * sqrt(sqrt(m22[iiP + ii] * m22[jjP + jj]) / m22[kkP + kk]) * temp_kk +
2262 sqrt(m22[kkP + kk] * sqrt(m22[jjP + jj] / (m22[iiP + ii] * m22[iiP + ii] * m22[iiP + ii]))) * temp_ii +
2263 sqrt(m22[kkP + kk] * sqrt(m22[iiP + ii] / (m22[jjP + jj] * m22[jjP + jj] * m22[jjP + jj]))) * temp_jj) / N;
2264 } else {
2265 // psi_ijkl
2266 double S5111 = 0.0;
2267 double S1511 = 0.0;
2268 double S1151 = 0.0;
2269 double S1115 = 0.0;
2270 double S2111 = 0.0;
2271 double S1211 = 0.0;
2272 double S1121 = 0.0;
2273 double S1112 = 0.0;
2274 double m1111 = 0.0;
2275 double m0111 = 0.0;
2276 double m1011 = 0.0;
2277 double m1101 = 0.0;
2278 double m1110 = 0.0;
2279 for (int tt = 0; tt < n; tt++) {
2280 S5111 += Xc2[iiN + tt] * Xc2[iiN + tt] * Xc[iiN + tt] * Xc[jjN + tt] * Xc[kkN + tt] * Xc[llN + tt];
2281 S1511 += Xc[iiN + tt] * Xc2[jjN + tt] * Xc2[jjN + tt] * Xc[jjN + tt] * Xc[kkN + tt] * Xc[llN + tt];
2282 S1151 += Xc[iiN + tt] * Xc[jjN + tt] * Xc2[kkN + tt] * Xc2[kkN + tt] * Xc[kkN + tt] * Xc[llN + tt];
2283 S1115 += Xc[iiN + tt] * Xc[jjN + tt] * Xc[kkN + tt] * Xc2[llN + tt] * Xc2[llN + tt] * Xc[llN + tt];
2284 S2111 += Xc2[iiN + tt] * Xc[jjN + tt] * Xc[kkN + tt] * Xc[llN + tt];
2285 S1211 += Xc[iiN + tt] * Xc2[jjN + tt] * Xc[kkN + tt] * Xc[llN + tt];
2286 S1121 += Xc[iiN + tt] * Xc[jjN + tt] * Xc2[kkN + tt] * Xc[llN + tt];
2287 S1112 += Xc[iiN + tt] * Xc[jjN + tt] * Xc[kkN + tt] * Xc2[llN + tt];
2288 m1111 += Xc[iiN + tt] * Xc[jjN + tt] * Xc[kkN + tt] * Xc[llN + tt];
2289 m0111 += Xc[jjN + tt] * Xc[kkN + tt] * Xc[llN + tt];
2290 m1011 += Xc[iiN + tt] * Xc[kkN + tt] * Xc[llN + tt];
2291 m1101 += Xc[iiN + tt] * Xc[jjN + tt] * Xc[llN + tt];
2292 m1110 += Xc[iiN + tt] * Xc[jjN + tt] * Xc[kkN + tt];
2293 }
2294 m1111 /= N;
2295 m0111 /= N;
2296 m1011 /= N;
2297 m1101 /= N;
2298 m1110 /= N;
2299
2300 double temp_ii = S5111 / N - m1111 / N * m22[iiP + ii] - 4.0 * S2111 / N * m21[iiP + ii] -
2301 m0111 * m41[iiP + ii] - m1011 * m41[jjP + ii] - m1101 * m41[kkP + ii] -
2302 m1110 * m41[llP + ii] + 4.0 * m1110 * m21[iiP + ii] * m11[llP + ii] +
2303 4.0 * m1101 * m21[iiP + ii] * m11[kkP + ii] + 4.0 * m1011 * m21[iiP + ii] * m11[jjP + ii] +
2304 4.0 * m0111 * m11[iiP + ii] * m21[iiP + ii];
2305 double temp_jj = S1511 / N - m1111 / N * m22[jjP + jj] - 4.0 * S1211 / N * m21[jjP + jj] -
2306 m0111 * m41[iiP + jj] - m1011 * m41[jjP + jj] - m1101 * m41[kkP + jj] -
2307 m1110 * m41[llP + jj] + 4.0 * m1110 * m21[jjP + jj] * m11[llP + jj] +
2308 4.0 * m1101 * m21[jjP + jj] * m11[kkP + jj] + 4.0 * m1011 * m21[jjP + jj] * m11[jjP + jj] +
2309 4.0 * m0111 * m11[iiP + jj] * m21[jjP + jj];
2310 double temp_kk = S1151 / N - m1111 / N * m22[kkP + kk] - 4.0 * S1121 / N * m21[kkP + kk] -
2311 m0111 * m41[iiP + kk] - m1011 * m41[jjP + kk] - m1101 * m41[kkP + kk] -
2312 m1110 * m41[llP + kk] + 4.0 * m1110 * m21[kkP + kk] * m11[llP + kk] +
2313 4.0 * m1101 * m21[kkP + kk] * m11[kkP + kk] + 4.0 * m1011 * m21[kkP + kk] * m11[jjP + kk] +
2314 4.0 * m0111 * m11[iiP + kk] * m21[kkP + kk];
2315 double temp_ll = S1115 / N - m1111 / N * m22[llP + ll] - 4.0 * S1112 / N * m21[llP + ll] -
2316 m0111 * m41[iiP + ll] - m1011 * m41[jjP + ll] - m1101 * m41[kkP + ll] -
2317 m1110 * m41[llP + ll] + 4.0 * m1110 * m21[llP + ll] * m11[llP + ll] +
2318 4.0 * m1101 * m21[llP + ll] * m11[kkP + ll] + 4.0 * m1011 * m21[llP + ll] * m11[jjP + ll] +
2319 4.0 * m0111 * m11[iiP + ll] * m21[llP + ll];
2320
2321 rCM4[0] += 6.0 * r7 * r5 *
2322 (sqrt(sqrt((m22[jjP + jj] * m22[kkP + kk] * m22[llP + ll]) / (m22[iiP + ii] * m22[iiP + ii] * m22[iiP + ii]))) * temp_ii +
2323 sqrt(sqrt((m22[iiP + ii] * m22[kkP + kk] * m22[llP + ll]) / (m22[jjP + jj] * m22[jjP + jj] * m22[jjP + jj]))) * temp_jj +
2324 sqrt(sqrt((m22[iiP + ii] * m22[jjP + jj] * m22[llP + ll]) / (m22[kkP + kk] * m22[kkP + kk] * m22[kkP + kk]))) * temp_kk +
2325 sqrt(sqrt((m22[iiP + ii] * m22[jjP + jj] * m22[kkP + kk]) / (m22[llP + ll] * m22[llP + ll] * m22[llP + ll]))) * temp_ll) / N;
2326 }
2327 }
2328 }
2329 } // loop ll
2330 } // loop kk
2331 } // loop jj
2332 } // loop ii
2333
2334 UNPROTECT(1);
2335 return CM4;
2336 }
2337
CM4_1F(SEXP XXc,SEXP XXc2,SEXP ffcobs,SEXP ffvar,SEXP ffskew,SEXP ffkurt,SEXP mm11,SEXP mm21,SEXP mm22,SEXP mm31,SEXP NN,SEXP PP)2338 SEXP CM4_1F(SEXP XXc, SEXP XXc2, SEXP ffcobs, SEXP ffvar, SEXP ffskew, SEXP ffkurt,
2339 SEXP mm11, SEXP mm21, SEXP mm22, SEXP mm31, SEXP NN, SEXP PP){
2340 /*
2341 arguments
2342 XXc : numeric vector with centered observations
2343 XXc2 : numeric vector with centered and squared observations
2344 ffcobs : centered factor observations
2345 ffvar : variance of the factor
2346 ffskew : third order central moment of the factor
2347 ffkurt : fourth order central moment of the factor
2348 mm11 : numeric vector of t(Xc) %*% Xc / NN
2349 mm21 : numeric vector of t(Xc^2) %*% Xc / NN
2350 mm22 : numeric vector of t(Xc^2) %*% Xc^2 / NN
2351 mm31 : numeric vector of t(Xc^3) %*% Xc / NN
2352 NN : integer, number of observations
2353 PP : integer, number of assets
2354
2355 Written by Dries Cornilly
2356 */
2357
2358 // // declare pointers for the vector
2359 double *Xc, *Xc2, *fcobs, *m11, *m21, *m22, *m31;
2360
2361 // use REAL() to access the C array inside the numeric vector passed in from R
2362 Xc = REAL(XXc);
2363 Xc2 = REAL(XXc2);
2364 fcobs = REAL(ffcobs);
2365 m11 = REAL(mm11);
2366 m21 = REAL(mm21);
2367 m22 = REAL(mm22);
2368 m31 = REAL(mm31);
2369 double N = asReal(NN);
2370 int n = asInteger(NN);
2371 int P = asInteger(PP);
2372 double fvar = asReal(ffvar);
2373 double fskew = asReal(ffskew);
2374 double fkurt = asReal(ffkurt);
2375
2376 // allocate and compute the covariance
2377 SEXP CM4 = PROTECT(allocVector(REALSXP, 1));
2378 double *rCM4 = REAL(CM4);
2379 rCM4[0] = 0.0;
2380
2381 // compute some helper variables for later on
2382 double fvar2 = fvar * fvar;
2383 double fvar3 = fvar2 * fvar;
2384 double fvar4 = fvar2 * fvar2;
2385 double fvar5 = fvar2 * fvar3;
2386 double covXf[P];
2387 double X1f2[P];
2388 double X1f4[P];
2389 double X11f1[P * P];
2390 for (int ii = 0; ii < P; ii++) {
2391 int iiN = ii * n;
2392 covXf[ii] = 0.0;
2393 X1f2[ii] = 0.0;
2394 X1f4[ii] = 0.0;
2395 for (int tt = 0; tt < n; tt++) {
2396 covXf[ii] += Xc[iiN + tt] * fcobs[tt];
2397 X1f2[ii] += Xc[iiN + tt] * fcobs[tt] * fcobs[tt];
2398 X1f4[ii] += Xc[iiN + tt] * fcobs[tt] * fcobs[tt] * fcobs[tt] * fcobs[tt];
2399 }
2400 covXf[ii] /= N;
2401 X1f2[ii] /= N;
2402 X1f4[ii] /= N;
2403 for (int jj = ii; jj < P; jj++) {
2404 int jjN = jj * n;
2405 double elem = 0.0;
2406 for (int tt = 0; tt < n; tt++) {
2407 elem += Xc[iiN + tt] * Xc[jjN + tt] * fcobs[tt];
2408 }
2409 elem /= N;
2410 X11f1[ii * P + jj] = elem;
2411 X11f1[jj * P + ii] = elem;
2412 }
2413 }
2414
2415 // loop over the unique indices i <= j <= k <= l
2416 for (int ii = 0; ii < P; ii++) {
2417 int iiP = ii * P;
2418 int iiN = ii * n;
2419 for (int jj = ii; jj < P; jj++) {
2420 int jjP = jj * P;
2421 int jjN = jj * n;
2422 for (int kk = jj; kk < P; kk++) {
2423 int kkP = kk * P;
2424 int kkN = kk * n;
2425 for (int ll = kk; ll < P; ll++) {
2426 int llP = ll * P;
2427 int llN = ll * n;
2428 if (ii == jj) {
2429 if (jj == kk) {
2430 if (kk == ll) {
2431 // psi_iiii
2432 double S8 = 0.0;
2433 double S5 = 0.0;
2434 for (int tt = 0; tt < n; tt++) {
2435 S8 += Xc2[iiN + tt] * Xc2[iiN + tt] * Xc2[iiN + tt] * Xc2[iiN + tt];
2436 S5 += Xc2[iiN + tt] * Xc2[iiN + tt] * Xc[iiN + tt];
2437 }
2438
2439 rCM4[0] += (S8 / N - m31[iiP + ii] * m31[iiP + ii] -
2440 8.0 * S5 / N * m21[iiP + ii] +
2441 16.0 * m11[iiP + ii] * m21[iiP + ii] * m21[iiP + ii]) / N;
2442 } else {
2443 // psi_iiil
2444 double S411 = 0.0;
2445 double S321 = 0.0;
2446 double S51 = 0.0;
2447 double S312 = 0.0;
2448 double S314 = 0.0;
2449 double S311 = 0.0;
2450 for (int tt = 0; tt < n; tt++) {
2451 S411 += Xc2[iiN + tt] * Xc2[iiN + tt] * Xc[llN + tt] * fcobs[tt];
2452 S321 += Xc2[iiN + tt] * Xc[iiN + tt] * Xc2[llN + tt] * fcobs[tt];
2453 S51 += Xc2[iiN + tt] * Xc2[iiN + tt] * Xc[iiN + tt] * Xc[llN + tt];
2454 S312 += Xc2[iiN + tt] * Xc[iiN + tt] * Xc[llN + tt] * fcobs[tt] * fcobs[tt];
2455 S314 += Xc2[iiN + tt] * Xc[iiN + tt] * Xc[llN + tt] * fcobs[tt] * fcobs[tt] * fcobs[tt] * fcobs[tt];
2456 S311 += Xc2[iiN + tt] * Xc[iiN + tt] * Xc[llN + tt] * fcobs[tt];
2457 }
2458
2459 double temp_i0 = S411 / N - m31[llP + ii] * covXf[ii] - 3.0 * m21[llP + ii] * X11f1[iiP + ii] -
2460 m21[iiP + ii] * X11f1[llP + ii];
2461 double temp_l0 = S321 / N - m31[llP + ii] * covXf[ll] - 3.0 * m21[llP + ii] * X11f1[llP + ii] -
2462 m21[iiP + ii] * X11f1[llP + ll];
2463 double temp_ii = S51 / N - m31[llP + ii] * m11[iiP + ii] - 4.0 * m21[llP + ii] * m21[iiP + ii];
2464 double temp_var = S312 / N - m31[llP + ii] * fvar - 3.0 * m21[llP + ii] * X1f2[ii] -
2465 m21[iiP + ii] * X1f2[ll];
2466 double temp_kurt = S314 / N - m31[llP + ii] * fkurt - 4.0 * S311 / N * fkurt +
2467 4.0 * fskew * (3.0 * m21[llP + ii] * covXf[ii] + m21[iiP + ii] * covXf[ll]);
2468
2469 rCM4[0] += 4.0 * ((3.0 * covXf[ii] * covXf[ii] * covXf[ll] * fkurt / fvar4 -
2470 9.0 * covXf[ii] * covXf[ii] * covXf[ll] / fvar2 + 3.0 * covXf[ll] * m11[iiP + ii] / fvar) * temp_i0 +
2471 (covXf[ii] * covXf[ii] * covXf[ii] * fkurt / fvar4 +
2472 3.0 * covXf[ii] * (m11[iiP + ii] - covXf[ii] * covXf[ii] / fvar) / fvar) * temp_l0 +
2473 3.0 * covXf[ii] * covXf[ll] / fvar * temp_ii +
2474 (-4.0 * covXf[ii] * covXf[ii] * covXf[ii] * covXf[ll] * fkurt / fvar5 -
2475 3.0 * covXf[ii] * covXf[ll] * m11[iiP + ii] / fvar2 +
2476 6.0 * covXf[ii] * covXf[ii] * covXf[ii] * covXf[ll] / fvar3) * temp_var +
2477 covXf[ii] * covXf[ii] * covXf[ii] * covXf[ll] * temp_kurt / fvar4) / N;
2478 }
2479 } else {
2480 if (kk == ll) {
2481 // psi_iikk
2482 double S321 = 0.0;
2483 double S231 = 0.0;
2484 double S420 = 0.0;
2485 double S240 = 0.0;
2486 double S222 = 0.0;
2487 double S224 = 0.0;
2488 for (int tt = 0; tt < n; tt++) {
2489 S321 += Xc2[iiN + tt] * Xc[iiN + tt] * Xc2[kkN + tt] * fcobs[tt];
2490 S231 += Xc2[iiN + tt] * Xc2[kkN + tt] * Xc[kkN + tt] * fcobs[tt];
2491 S420 += Xc2[iiN + tt] * Xc2[iiN + tt] * Xc2[kkN + tt];
2492 S240 += Xc2[iiN + tt] * Xc2[kkN + tt] * Xc2[kkN + tt];
2493 S222 += Xc2[iiN + tt] * Xc2[kkN + tt] * fcobs[tt] * fcobs[tt];
2494 S224 += Xc2[iiN + tt] * Xc2[kkN + tt] * fcobs[tt] * fcobs[tt] * fcobs[tt] * fcobs[tt];
2495 }
2496
2497 double temp_i0 = S321 / N - m22[kkP + ii] * covXf[ii] - 2.0 * m21[iiP + kk] * X11f1[iiP + ii] -
2498 2.0 * m21[kkP + ii] * X11f1[kkP + ii];
2499 double temp_k0 = S231 / N - m22[iiP + kk] * covXf[kk] - 2.0 * m21[kkP + ii] * X11f1[kkP + kk] -
2500 2.0 * m21[iiP + kk] * X11f1[iiP + kk];
2501 double temp_ii = S420 / N - m22[kkP + ii] * m11[iiP + ii] - 2.0 * m21[iiP + kk] * m21[iiP + ii] -
2502 2.0 * m21[kkP + ii] * m21[kkP + ii];
2503 double temp_kk = S240 / N - m22[iiP + kk] * m11[kkP + kk] - 2.0 * m21[kkP + ii] * m21[kkP + kk] -
2504 2.0 * m21[iiP + kk] * m21[iiP + kk];
2505 double temp_var = S222 / N - m22[kkP + ii] * fvar - 2.0 * m21[iiP + kk] * X1f2[ii] -
2506 2.0 * m21[kkP + ii] * X1f2[kk];
2507 double temp_kurt = S224 / N - m22[kkP + ii] * fkurt - 2.0 * m21[iiP + kk] * X1f4[ii] -
2508 2.0 * m21[kkP + ii] * X1f4[kk] + 8.0 * fskew * (m21[iiP + kk] * covXf[ii] + m21[kkP + ii] * covXf[kk]);
2509
2510 rCM4[0] += 6.0 * ((2.0 * covXf[ii] * covXf[kk] * covXf[kk] * fkurt / fvar4 -
2511 2.0 * covXf[ii] * covXf[kk] * covXf[kk] / fvar2) * temp_i0 +
2512 (2.0 * covXf[ii] * covXf[ii] * covXf[kk] * fkurt / fvar4 -
2513 2.0 * covXf[ii] * covXf[ii] * covXf[kk] / fvar2) * temp_k0 +
2514 m11[kkP + kk] * temp_ii + m11[iiP + ii] * temp_kk +
2515 (-4.0 * covXf[ii] * covXf[ii] * covXf[kk] * covXf[kk] * fkurt / fvar5 +
2516 2.0 * covXf[ii] * covXf[ii] * covXf[kk] * covXf[kk] / fvar3) * temp_var +
2517 covXf[ii] * covXf[ii] * covXf[kk] * covXf[kk] * temp_kurt / fvar4) / N;
2518 } else {
2519 // psi_iikl
2520 double S3111 = 0.0;
2521 double S2211 = 0.0;
2522 double S2121 = 0.0;
2523 double S4110 = 0.0;
2524 double S2112 = 0.0;
2525 double S2114 = 0.0;
2526 double S2111 = 0.0;
2527 double m2110 = 0.0;
2528 double m1110 = 0.0;
2529 for (int tt = 0; tt < n; tt++) {
2530 S3111 += Xc2[iiN + tt] * Xc[iiN + tt] * Xc[kkN + tt] * Xc[llN + tt] * fcobs[tt];
2531 S2211 += Xc2[iiN + tt] * Xc2[kkN + tt] * Xc[llN + tt] * fcobs[tt];
2532 S2121 += Xc2[iiN + tt] * Xc[kkN + tt] * Xc2[llN + tt] * fcobs[tt];
2533 S4110 += Xc2[iiN + tt] * Xc2[iiN + tt] * Xc[kkN + tt] * Xc[llN + tt];
2534 S2112 += Xc2[iiN + tt] * Xc[kkN + tt] * Xc[llN + tt] * fcobs[tt] * fcobs[tt];
2535 S2114 += Xc2[iiN + tt] * Xc[kkN + tt] * Xc[llN + tt] * fcobs[tt] * fcobs[tt] * fcobs[tt] * fcobs[tt];
2536 S2111 += Xc2[iiN + tt] * Xc[kkN + tt] * Xc[llN + tt] * fcobs[tt];
2537 m2110 += Xc2[iiN + tt] * Xc[kkN + tt] * Xc[llN + tt];
2538 m1110 += Xc[iiN + tt] * Xc[kkN + tt] * Xc[llN + tt];
2539 }
2540 m2110 /= N;
2541 m1110 /= N;
2542
2543 double temp_i0 = S3111 / N - m2110 * covXf[ii] - 2.0 * m1110 * X11f1[iiP + ii] -
2544 m21[llP + ii] * X11f1[kkP + ii] - m21[kkP + ii] * X11f1[llP + ii];
2545 double temp_k0 = S2211 / N - m2110 * covXf[kk] - 2.0 * m1110 * X11f1[kkP + ii] -
2546 m21[llP + ii] * X11f1[kkP + kk] - m21[kkP + ii] * X11f1[llP + kk];
2547 double temp_l0 = S2121 / N - m2110 * covXf[ll] - 2.0 * m1110 * X11f1[llP + ii] -
2548 m21[kkP + ii] * X11f1[llP + ll] - m21[llP + ii] * X11f1[kkP + ll];
2549 double temp_ii = S4110 / N - m2110 * m11[iiP + ii] - 2.0 * m1110 * m21[iiP + ii] -
2550 2.0 * m21[llP + ii] * m21[kkP + ii];
2551 double temp_var = S2112 / N - m2110 * fvar - 2.0 * m1110 * X1f2[ii] -
2552 m21[llP + ii] * X1f2[kk] - m21[kkP + ii] * X1f2[ll];
2553 double temp_kurt = S2114 / N - m2110 * fkurt - 2.0 * m1110 * X1f4[ii] -
2554 m21[llP + ii] * X1f4[kk] - m21[kkP + ii] * X1f4[ll] - 4.0 * S2111 / N * fskew +
2555 4.0 * fskew * (2.0 * m1110 * covXf[ii] + m21[llP + ii] * covXf[kk] + m21[kkP + ii] * covXf[ll]);
2556
2557 rCM4[0] += 12.0 * ((2.0 * covXf[ii] * covXf[kk] * covXf[ll] * fkurt / fvar4 -
2558 2.0 * covXf[ii] * covXf[kk] * covXf[ll] / fvar2) * temp_i0 +
2559 (covXf[ii] * covXf[ii] * covXf[ll] * fkurt / fvar4 +
2560 covXf[ll] * (m11[iiP + ii] - covXf[ii] * covXf[ii] / fvar) / fvar) * temp_k0 +
2561 (covXf[ii] * covXf[ii] * covXf[kk] * fkurt / fvar4 +
2562 covXf[kk] * (m11[iiP + ii] - covXf[ii] * covXf[ii] / fvar) / fvar) * temp_l0 +
2563 covXf[kk] * covXf[ll] / fvar * temp_ii +
2564 (-4.0 * covXf[ii] * covXf[ii] * covXf[kk] * covXf[ll] * fkurt / fvar5 +
2565 2.0 * covXf[kk] * covXf[ll] * covXf[ii] * covXf[ii] / fvar3 -
2566 covXf[kk] * covXf[ll] * m11[iiP + ii] / fvar2) * temp_var +
2567 covXf[ii] * covXf[ii] * covXf[kk] * covXf[ll] * temp_kurt / fvar4) / N;
2568 }
2569 }
2570 } else {
2571 if (jj == kk) {
2572 if (kk == ll) {
2573 // psi_ijjj
2574 double S411 = 0.0;
2575 double S321 = 0.0;
2576 double S51 = 0.0;
2577 double S312 = 0.0;
2578 double S314 = 0.0;
2579 double S311 = 0.0;
2580 for (int tt = 0; tt < n; tt++) {
2581 S411 += Xc2[jjN + tt] * Xc2[jjN + tt] * Xc[iiN + tt] * fcobs[tt];
2582 S321 += Xc2[jjN + tt] * Xc[jjN + tt] * Xc2[iiN + tt] * fcobs[tt];
2583 S51 += Xc2[jjN + tt] * Xc2[jjN + tt] * Xc[jjN + tt] * Xc[iiN + tt];
2584 S312 += Xc2[jjN + tt] * Xc[jjN + tt] * Xc[iiN + tt] * fcobs[tt] * fcobs[tt];
2585 S314 += Xc2[jjN + tt] * Xc[jjN + tt] * Xc[iiN + tt] * fcobs[tt] * fcobs[tt] * fcobs[tt] * fcobs[tt];
2586 S311 += Xc2[jjN + tt] * Xc[jjN + tt] * Xc[iiN + tt] * fcobs[tt];
2587 }
2588
2589 double temp_j0 = S411 / N - m31[iiP + jj] * covXf[jj] - 3.0 * m21[iiP + jj] * X11f1[jjP + jj] -
2590 m21[jjP + jj] * X11f1[iiP + jj];
2591 double temp_i0 = S321 / N - m31[iiP + jj] * covXf[ii] - 3.0 * m21[iiP + jj] * X11f1[iiP + jj] -
2592 m21[jjP + jj] * X11f1[iiP + ii];
2593 double temp_jj = S51 / N - m31[iiP + jj] * m11[jjP + jj] - 4.0 * m21[iiP + jj] * m21[jjP + jj];
2594 double temp_var = S312 / N - m31[iiP + jj] * fvar - 3.0 * m21[iiP + jj] * X1f2[jj] -
2595 m21[jjP + jj] * X1f2[ii];
2596 double temp_kurt = S314 / N - m31[iiP + jj] * fkurt - 4.0 * S311 / N * fkurt +
2597 4.0 * fskew * (3.0 * m21[iiP + jj] * covXf[jj] + m21[jjP + jj] * covXf[ii]);
2598
2599 rCM4[0] += 4.0 * ((3.0 * covXf[jj] * covXf[jj] * covXf[ii] * fkurt / fvar4 -
2600 9.0 * covXf[jj] * covXf[jj] * covXf[ii] / fvar2 + 3.0 * covXf[ii] * m11[jjP + jj] / fvar) * temp_j0 +
2601 (covXf[jj] * covXf[jj] * covXf[jj] * fkurt / fvar4 +
2602 3.0 * covXf[jj] * (m11[jjP + jj] - covXf[jj] * covXf[jj] / fvar) / fvar) * temp_i0 +
2603 3.0 * covXf[jj] * covXf[ii] / fvar * temp_jj +
2604 (-4.0 * covXf[jj] * covXf[jj] * covXf[jj] * covXf[ii] * fkurt / fvar5 -
2605 3.0 * covXf[jj] * covXf[ii] * m11[jjP + jj] / fvar2 +
2606 6.0 * covXf[jj] * covXf[jj] * covXf[jj] * covXf[ii] / fvar3) * temp_var +
2607 covXf[jj] * covXf[jj] * covXf[jj] * covXf[ii] * temp_kurt / fvar4) / N;
2608 } else {
2609 // psi_ijjl
2610 double S3111 = 0.0;
2611 double S2211 = 0.0;
2612 double S2121 = 0.0;
2613 double S4110 = 0.0;
2614 double S2112 = 0.0;
2615 double S2114 = 0.0;
2616 double S2111 = 0.0;
2617 double m2110 = 0.0;
2618 double m1110 = 0.0;
2619 for (int tt = 0; tt < n; tt++) {
2620 S3111 += Xc2[jjN + tt] * Xc[jjN + tt] * Xc[iiN + tt] * Xc[llN + tt] * fcobs[tt];
2621 S2211 += Xc2[jjN + tt] * Xc2[iiN + tt] * Xc[llN + tt] * fcobs[tt];
2622 S2121 += Xc2[jjN + tt] * Xc[iiN + tt] * Xc2[llN + tt] * fcobs[tt];
2623 S4110 += Xc2[jjN + tt] * Xc2[jjN + tt] * Xc[iiN + tt] * Xc[llN + tt];
2624 S2112 += Xc2[jjN + tt] * Xc[iiN + tt] * Xc[llN + tt] * fcobs[tt] * fcobs[tt];
2625 S2114 += Xc2[jjN + tt] * Xc[iiN + tt] * Xc[llN + tt] * fcobs[tt] * fcobs[tt] * fcobs[tt] * fcobs[tt];
2626 S2111 += Xc2[jjN + tt] * Xc[iiN + tt] * Xc[llN + tt] * fcobs[tt];
2627 m2110 += Xc2[jjN + tt] * Xc[iiN + tt] * Xc[llN + tt];
2628 m1110 += Xc[jjN + tt] * Xc[iiN + tt] * Xc[llN + tt];
2629 }
2630 m2110 /= N;
2631 m1110 /= N;
2632
2633 double temp_j0 = S3111 / N - m2110 * covXf[jj] - 2.0 * m1110 * X11f1[jjP + jj] -
2634 m21[llP + jj] * X11f1[iiP + jj] - m21[iiP + jj] * X11f1[llP + jj];
2635 double temp_i0 = S2211 / N - m2110 * covXf[ii] - 2.0 * m1110 * X11f1[iiP + jj] -
2636 m21[llP + jj] * X11f1[iiP + ii] - m21[iiP + jj] * X11f1[llP + ii];
2637 double temp_l0 = S2121 / N - m2110 * covXf[ll] - 2.0 * m1110 * X11f1[llP + jj] -
2638 m21[iiP + jj] * X11f1[llP + ll] - m21[llP + jj] * X11f1[iiP + ll];
2639 double temp_jj = S4110 / N - m2110 * m11[jjP + jj] - 2.0 * m1110 * m21[jjP + jj] -
2640 2.0 * m21[llP + jj] * m21[iiP + jj];
2641 double temp_var = S2112 / N - m2110 * fvar - 2.0 * m1110 * X1f2[jj] -
2642 m21[llP + jj] * X1f2[ii] - m21[iiP + jj] * X1f2[ll];
2643 double temp_kurt = S2114 / N - m2110 * fkurt - 2.0 * m1110 * X1f4[jj] -
2644 m21[llP + jj] * X1f4[ii] - m21[iiP + jj] * X1f4[ll] - 4.0 * S2111 / N * fskew +
2645 4.0 * fskew * (2.0 * m1110 * covXf[jj] + m21[llP + jj] * covXf[ii] + m21[iiP + jj] * covXf[ll]);
2646
2647 rCM4[0] += 12.0 * ((2.0 * covXf[jj] * covXf[ii] * covXf[ll] * fkurt / fvar4 -
2648 2.0 * covXf[jj] * covXf[ii] * covXf[ll] / fvar2) * temp_j0 +
2649 (covXf[jj] * covXf[jj] * covXf[ll] * fkurt / fvar4 +
2650 covXf[ll] * (m11[jjP + jj] - covXf[jj] * covXf[jj] / fvar) / fvar) * temp_i0 +
2651 (covXf[jj] * covXf[jj] * covXf[ii] * fkurt / fvar4 +
2652 covXf[ii] * (m11[jjP + jj] - covXf[jj] * covXf[jj] / fvar) / fvar) * temp_l0 +
2653 covXf[ii] * covXf[ll] / fvar * temp_jj +
2654 (-4.0 * covXf[jj] * covXf[jj] * covXf[ii] * covXf[ll] * fkurt / fvar5 +
2655 2.0 * covXf[ii] * covXf[ll] * covXf[jj] * covXf[jj] / fvar3 -
2656 covXf[ii] * covXf[ll] * m11[jjP + jj] / fvar2) * temp_var +
2657 covXf[jj] * covXf[jj] * covXf[ii] * covXf[ll] * temp_kurt / fvar4) / N;
2658 }
2659 } else {
2660 if (kk == ll) {
2661 // psi_ijkk
2662 double S3111 = 0.0;
2663 double S2211 = 0.0;
2664 double S2121 = 0.0;
2665 double S4110 = 0.0;
2666 double S2112 = 0.0;
2667 double S2114 = 0.0;
2668 double S2111 = 0.0;
2669 double m2110 = 0.0;
2670 double m1110 = 0.0;
2671 for (int tt = 0; tt < n; tt++) {
2672 S3111 += Xc2[kkN + tt] * Xc[kkN + tt] * Xc[iiN + tt] * Xc[jjN + tt] * fcobs[tt];
2673 S2211 += Xc2[kkN + tt] * Xc2[iiN + tt] * Xc[jjN + tt] * fcobs[tt];
2674 S2121 += Xc2[kkN + tt] * Xc[iiN + tt] * Xc2[jjN + tt] * fcobs[tt];
2675 S4110 += Xc2[kkN + tt] * Xc2[kkN + tt] * Xc[iiN + tt] * Xc[jjN + tt];
2676 S2112 += Xc2[kkN + tt] * Xc[iiN + tt] * Xc[jjN + tt] * fcobs[tt] * fcobs[tt];
2677 S2114 += Xc2[kkN + tt] * Xc[iiN + tt] * Xc[jjN + tt] * fcobs[tt] * fcobs[tt] * fcobs[tt] * fcobs[tt];
2678 S2111 += Xc2[kkN + tt] * Xc[iiN + tt] * Xc[jjN + tt] * fcobs[tt];
2679 m2110 += Xc2[kkN + tt] * Xc[iiN + tt] * Xc[jjN + tt];
2680 m1110 += Xc[kkN + tt] * Xc[iiN + tt] * Xc[jjN + tt];
2681 }
2682 m2110 /= N;
2683 m1110 /= N;
2684
2685 double temp_k0 = S3111 / N - m2110 * covXf[kk] - 2.0 * m1110 * X11f1[kkP + kk] -
2686 m21[jjP + kk] * X11f1[iiP + kk] - m21[iiP + kk] * X11f1[jjP + kk];
2687 double temp_i0 = S2211 / N - m2110 * covXf[ii] - 2.0 * m1110 * X11f1[iiP + kk] -
2688 m21[jjP + kk] * X11f1[iiP + ii] - m21[iiP + kk] * X11f1[jjP + ii];
2689 double temp_j0 = S2121 / N - m2110 * covXf[jj] - 2.0 * m1110 * X11f1[jjP + kk] -
2690 m21[iiP + kk] * X11f1[jjP + jj] - m21[jjP + kk] * X11f1[iiP + jj];
2691 double temp_kk = S4110 / N - m2110 * m11[kkP + kk] - 2.0 * m1110 * m21[kkP + kk] -
2692 2.0 * m21[jjP + kk] * m21[iiP + kk];
2693 double temp_var = S2112 / N - m2110 * fvar - 2.0 * m1110 * X1f2[kk] -
2694 m21[jjP + kk] * X1f2[ii] - m21[iiP + kk] * X1f2[jj];
2695 double temp_kurt = S2114 / N - m2110 * fkurt - 2.0 * m1110 * X1f4[kk] -
2696 m21[jjP + kk] * X1f4[ii] - m21[iiP + kk] * X1f4[jj] - 4.0 * S2111 / N * fskew +
2697 4.0 * fskew * (2.0 * m1110 * covXf[kk] + m21[jjP + kk] * covXf[ii] + m21[iiP + kk] * covXf[jj]);
2698
2699 rCM4[0] += 12.0 * ((2.0 * covXf[kk] * covXf[ii] * covXf[jj] * fkurt / fvar4 -
2700 2.0 * covXf[kk] * covXf[ii] * covXf[jj] / fvar2) * temp_k0 +
2701 (covXf[kk] * covXf[kk] * covXf[jj] * fkurt / fvar4 +
2702 covXf[jj] * (m11[kkP + kk] - covXf[kk] * covXf[kk] / fvar) / fvar) * temp_i0 +
2703 (covXf[kk] * covXf[kk] * covXf[ii] * fkurt / fvar4 +
2704 covXf[ii] * (m11[kkP + kk] - covXf[kk] * covXf[kk] / fvar) / fvar) * temp_j0 +
2705 covXf[ii] * covXf[jj] / fvar * temp_kk +
2706 (-4.0 * covXf[kk] * covXf[kk] * covXf[ii] * covXf[jj] * fkurt / fvar5 +
2707 2.0 * covXf[ii] * covXf[jj] * covXf[kk] * covXf[kk] / fvar3 -
2708 covXf[ii] * covXf[jj] * m11[kkP + kk] / fvar2) * temp_var +
2709 covXf[kk] * covXf[kk] * covXf[ii] * covXf[jj] * temp_kurt / fvar4) / N;
2710 } else {
2711 // psi_ijkl
2712 double S11114 = 0.0;
2713 double S11112 = 0.0;
2714 double S11111 = 0.0;
2715 double S21111 = 0.0;
2716 double S12111 = 0.0;
2717 double S11211 = 0.0;
2718 double S11121 = 0.0;
2719 double m1111 = 0.0;
2720 double m01110 = 0.0;
2721 double m10110 = 0.0;
2722 double m11010 = 0.0;
2723 double m11100 = 0.0;
2724 for (int tt = 0; tt < n; tt++) {
2725 S11114 += Xc[iiN + tt] * Xc[jjN + tt] * Xc[kkN + tt] * Xc[llN + tt] *
2726 fcobs[tt] * fcobs[tt] * fcobs[tt] * fcobs[tt];
2727 S11112 += Xc[iiN + tt] * Xc[jjN + tt] * Xc[kkN + tt] * Xc[llN + tt] * fcobs[tt] * fcobs[tt];
2728 S11111 += Xc[iiN + tt] * Xc[jjN + tt] * Xc[kkN + tt] * Xc[llN + tt] * fcobs[tt];
2729 S21111 += Xc2[iiN + tt] * Xc[jjN + tt] * Xc[kkN + tt] * Xc[llN + tt] * fcobs[tt];
2730 S12111 += Xc[iiN + tt] * Xc2[jjN + tt] * Xc[kkN + tt] * Xc[llN + tt] * fcobs[tt];
2731 S11211 += Xc[iiN + tt] * Xc[jjN + tt] * Xc2[kkN + tt] * Xc[llN + tt] * fcobs[tt];
2732 S11121 += Xc[iiN + tt] * Xc[jjN + tt] * Xc[kkN + tt] * Xc2[llN + tt] * fcobs[tt];
2733 m1111 += Xc[iiN + tt] * Xc[jjN + tt] * Xc[kkN + tt] * Xc[llN + tt];
2734 m01110 += Xc[jjN + tt] * Xc[kkN + tt] * Xc[llN + tt];
2735 m10110 += Xc[iiN + tt] * Xc[kkN + tt] * Xc[llN + tt];
2736 m11010 += Xc[iiN + tt] * Xc[jjN + tt] * Xc[llN + tt];
2737 m11100 += Xc[iiN + tt] * Xc[jjN + tt] * Xc[kkN + tt];
2738 }
2739 m1111 /= N;
2740 m01110 /= N;
2741 m10110 /= N;
2742 m11010 /= N;
2743 m11100 /= N;
2744
2745 double temp_ii = S21111 / N - m1111 * covXf[ii] - m01110 * X11f1[iiP + ii] - m10110 * X11f1[jjP + ii] -
2746 m11010 * X11f1[kkP + ii] - m11100 * X11f1[llP + ii];
2747 double temp_jj = S12111 / N - m1111 * covXf[jj] - m01110 * X11f1[iiP + jj] - m10110 * X11f1[jjP + jj] -
2748 m11010 * X11f1[kkP + jj] - m11100 * X11f1[llP + jj];
2749 double temp_kk = S11211 / N - m1111 * covXf[kk] - m01110 * X11f1[iiP + kk] - m10110 * X11f1[jjP + kk] -
2750 m11010 * X11f1[kkP + kk] - m11100 * X11f1[llP + kk];
2751 double temp_ll = S11121 / N - m1111 * covXf[ll] - m01110 * X11f1[iiP + ll] - m10110 * X11f1[jjP + ll] -
2752 m11010 * X11f1[kkP + ll] - m11100 * X11f1[llP + ll];
2753 double temp_kurt = S11114 / N - m1111 * fkurt - m01110 * X1f4[ii] - m10110 * X1f4[jj] -
2754 m11010 * X1f4[kk] - m11100 * X1f4[ll] - 4.0 * S11111 / N * fskew +
2755 4.0 * fskew * (m01110 * covXf[ii] + m10110 * covXf[jj] + m11010 * covXf[kk] + m11100 * covXf[ll]);
2756 double temp_var = S11112 / N - m1111 * fvar - m01110 * X1f2[ii] - m10110 * X1f2[jj] -
2757 m11010 * X1f2[kk] - m11100 * X1f2[ll];
2758
2759 rCM4[0] += 24.0 * ((covXf[jj] * covXf[kk] * covXf[ll] * temp_ii +
2760 covXf[ii] * covXf[kk] * covXf[ll] * temp_jj +
2761 covXf[ii] * covXf[jj] * covXf[ll] * temp_kk +
2762 covXf[ii] * covXf[jj] * covXf[kk] * temp_ll) * fkurt +
2763 covXf[ii] * covXf[jj] * covXf[kk] * covXf[ll] * temp_kurt -
2764 4.0 * covXf[ii] * covXf[jj] * covXf[kk] * covXf[ll] * fkurt * temp_var / fvar) / (N * fvar4);
2765 }
2766 }
2767 }
2768 } // loop ll
2769 } // loop kk
2770 } // loop jj
2771 } // loop ii
2772
2773 UNPROTECT(1);
2774 return CM4;
2775 }
2776