1/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3/*
4 Copyright (C) 2003 Ferdinando Ametrano
5 Copyright (C) 2003 RiskMap srl
6
7 This file is part of QuantLib, a free-software/open-source library
8 for financial quantitative analysts and developers - http://quantlib.org/
9
10 QuantLib is free software: you can redistribute it and/or modify it
11 under the terms of the QuantLib license. You should have received a
12 copy of the license along with this program; if not, please email
13 <quantlib-dev@lists.sf.net>. The license is also available online at
14 <http://quantlib.org/license.shtml>.
15
16 This program is distributed in the hope that it will be useful, but WITHOUT
17 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
18 FOR A PARTICULAR PURPOSE. See the license for more details.
19*/
20
21#include <ql/math/statistics/generalstatistics.hpp>
22
23namespace QuantLib {
24
25 Real GeneralStatistics::weightSum() const {
26 Real result = 0.0;
27 std::vector<std::pair<Real,Real> >::const_iterator it;
28 for (it=samples_.begin(); it!=samples_.end(); ++it) {
29 result += it->second;
30 }
31 return result;
32 }
33
34 Real GeneralStatistics::mean() const {
35 Size N = samples();
36 QL_REQUIRE(N != 0, "empty sample set");
37 return expectationValue(f: [](Real x) { return x; }).first;
38 }
39
40 Real GeneralStatistics::variance() const {
41 Size N = samples();
42 QL_REQUIRE(N > 1,
43 "sample number <=1, unsufficient");
44 // Subtract the mean and square. Repeat on the whole range.
45 // Hopefully, the whole thing will be inlined in a single loop.
46 Real m = mean();
47 Real s2 = expectationValue(f: [=](Real x) -> Real {
48 Real d = x - m;
49 return d * d;
50 }).first;
51 return s2*N/(N-1.0);
52 }
53
54 Real GeneralStatistics::skewness() const {
55 Size N = samples();
56 QL_REQUIRE(N > 2,
57 "sample number <=2, unsufficient");
58
59 Real m = mean();
60 Real X = expectationValue(f: [=](Real x) -> Real {
61 Real d = x - m;
62 return d * d * d;
63 }).first;
64 Real sigma = standardDeviation();
65
66 return (X/(sigma*sigma*sigma))*(N/(N-1.0))*(N/(N-2.0));
67 }
68
69 Real GeneralStatistics::kurtosis() const {
70 Size N = samples();
71 QL_REQUIRE(N > 3,
72 "sample number <=3, unsufficient");
73
74 Real m = mean();
75 Real X = expectationValue(f: [=](Real x) -> Real {
76 Real d = x - m;
77 Real d2 = d * d;
78 return d2 * d2;
79 }).first;
80 Real sigma2 = variance();
81
82 Real c1 = (N/(N-1.0)) * (N/(N-2.0)) * ((N+1.0)/(N-3.0));
83 Real c2 = 3.0 * ((N-1.0)/(N-2.0)) * ((N-1.0)/(N-3.0));
84
85 return c1*(X/(sigma2*sigma2))-c2;
86 }
87
88 Real GeneralStatistics::percentile(Real percent) const {
89
90 QL_REQUIRE(percent > 0.0 && percent <= 1.0,
91 "percentile (" << percent << ") must be in (0.0, 1.0]");
92
93 Real sampleWeight = weightSum();
94 QL_REQUIRE(sampleWeight>0.0,
95 "empty sample set");
96
97 sort();
98
99 std::vector<std::pair<Real,Real> >::iterator k, l;
100 k = samples_.begin();
101 l = samples_.end()-1;
102 /* the sum of weight is non null, therefore there's
103 at least one sample */
104 Real integral = k->second, target = percent*sampleWeight;
105 while (integral < target && k != l) {
106 ++k;
107 integral += k->second;
108 }
109 return k->first;
110 }
111
112 Real GeneralStatistics::topPercentile(Real percent) const {
113
114 QL_REQUIRE(percent > 0.0 && percent <= 1.0,
115 "percentile (" << percent << ") must be in (0.0, 1.0]");
116
117 Real sampleWeight = weightSum();
118 QL_REQUIRE(sampleWeight > 0.0,
119 "empty sample set");
120
121 sort();
122
123 std::vector<std::pair<Real,Real> >::reverse_iterator k, l;
124 k = samples_.rbegin();
125 l = samples_.rend()-1;
126 /* the sum of weight is non null, therefore there's
127 at least one sample */
128 Real integral = k->second, target = percent*sampleWeight;
129 while (integral < target && k != l) {
130 ++k;
131 integral += k->second;
132 }
133 return k->first;
134 }
135
136}
137

source code of quantlib/ql/math/statistics/generalstatistics.cpp