31 const Real tol_level = 1.0e-8;
33 class Garch11Constraint :
public Constraint {
35 class Impl final :
public Constraint::Impl {
38 Impl (
Real gammaLower,
Real gammaUpper)
40 bool test(
const Array& x)
const override {
41 QL_REQUIRE(x.size() >= 3,
"size of parameters vector < 3");
42 return x[0] > 0 && x[1] >= 0 && x[2] >= 0
48 Garch11Constraint(
Real gammaLower,
Real gammaUpper)
49 : Constraint(ext::shared_ptr<Constraint::Impl>(
50 new Garch11Constraint::Impl(gammaLower, gammaUpper))) {}
54 class Garch11CostFunction :
public CostFunction {
56 explicit Garch11CostFunction (
const std::vector<Volatility> &);
57 Real value(
const Array& x)
const override;
58 Array values(
const Array& x)
const override;
59 void gradient(Array& grad,
const Array& x)
const override;
60 Real valueAndGradient(Array& grad,
const Array& x)
const override;
63 const std::vector<Volatility> &
r2_;
66 Garch11CostFunction::Garch11CostFunction(
67 const std::vector<Volatility> &r2)
70 Real Garch11CostFunction::value(
const Array& x)
const {
75 sigma2 = x[0] + x[1] * u2 + x[2] * sigma2;
77 retval += std::log(sigma2) + u2 / sigma2;
79 return retval / (2.0*
r2_.size());
82 Array Garch11CostFunction::values(
const Array& x)
const {
83 Array retval (
r2_.size());
88 sigma2 = x[0] + x[1] * u2 + x[2] * sigma2;
90 retval[i++] = (std::log(sigma2) + u2 / sigma2)/(2.0*
r2_.size());
95 void Garch11CostFunction::gradient(Array& grad,
const Array& x)
const {
96 std::fill (grad.begin(), grad.end(), 0.0);
99 Real sigma2prev = sigma2;
102 for (
auto r2 :
r2_) {
103 sigma2 = x[0] + x[1] * u2 + x[2] * sigma2;
105 Real w = (sigma2 - u2) / (sigma2*sigma2);
107 grad[1] += u2prev * w;
108 grad[2] += sigma2prev * w;
112 std::transform(grad.begin(), grad.end(), grad.begin(),
113 [=](Real x) -> Real { return x / norm; });
116 Real Garch11CostFunction::valueAndGradient(Array& grad,
117 const Array& x)
const {
118 std::fill (grad.begin(), grad.end(), 0.0);
122 Real sigma2prev = sigma2;
125 for (
auto r2 :
r2_) {
126 sigma2 = x[0] + x[1] * u2 + x[2] * sigma2;
128 retval += std::log(sigma2) + u2 / sigma2;
129 Real w = (sigma2 - u2) / (sigma2*sigma2);
131 grad[1] += u2prev * w;
132 grad[2] += sigma2prev * w;
136 std::transform(grad.begin(), grad.end(), grad.begin(),
137 [=](Real x) -> Real { return x / norm; });
138 return retval / norm;
142 class FitAcfProblem :
public LeastSquareProblem {
144 FitAcfProblem(Real A2, Array acf, std::vector<std::size_t> idx);
145 Size size()
override;
146 void targetAndValue(
const Array& x, Array& target, Array& fct2fit)
override;
147 void targetValueAndGradient(
const Array& x,
148 Matrix& grad_fct2fit,
150 Array& fct2fit)
override;
158 FitAcfProblem::FitAcfProblem(Real A2, Array acf, std::vector<std::size_t> idx)
161 Size FitAcfProblem::size() {
return idx_.size(); }
163 void FitAcfProblem::targetAndValue(
const Array& x, Array& target,
171 / (3*(1 - gamma*gamma));
172 target[1] =
acf_[1] / A4;
173 fct2fit[1] = gamma * (1 - fct2fit[0]) -
beta;
174 for (std::size_t i = 2; i <
idx_.size(); ++i) {
176 fct2fit[i] = std::pow(gamma, (
int)
idx_[i]-1)* fct2fit[1];
180 void FitAcfProblem::targetValueAndGradient(
const Array& x,
181 Matrix& grad_fct2fit,
189 Real w2 = (1 - gamma*gamma);
190 fct2fit[0] = w1 / (3*w2);
191 grad_fct2fit[0][0] = (2.0/3.0) * ((2*
beta-3*gamma)*w2 + 2*w1*gamma) / (w2*w2);
192 grad_fct2fit[0][1] = (4.0/3.0) * (gamma -
beta) / w2;
193 target[1] =
acf_[1] / A4;
194 fct2fit[1] = gamma * (1 - fct2fit[0]) -
beta;
195 grad_fct2fit[1][0] = (1 - fct2fit[0]) - gamma * grad_fct2fit[0][0];
196 grad_fct2fit[1][1] = -gamma * grad_fct2fit[0][1] - 1;
197 for (std::size_t i = 2; i <
idx_.size(); ++i) {
199 w1 = std::pow(gamma, (
int)
idx_[i]-1);
200 fct2fit[i] = w1 * fct2fit[1];
201 grad_fct2fit[i][0] = (
idx_[i]-1) * (w1/gamma)*fct2fit[1] + w1*grad_fct2fit[1][0];
202 grad_fct2fit[i][1] = w1 * grad_fct2fit[1][1];
207 class FitAcfConstraint :
public Constraint {
209 class Impl final :
public Constraint::Impl {
212 Impl(Real gammaLower, Real gammaUpper)
214 bool test(
const Array& x)
const override {
215 QL_REQUIRE(x.size() >= 2,
"size of parameters vector < 2");
217 && x[1] >= 0 && x[1] <= x[0];
221 FitAcfConstraint(Real gammaLower, Real gammaUpper)
222 : Constraint(ext::shared_ptr<Constraint::Impl>(
223 new FitAcfConstraint::Impl(gammaLower, gammaUpper))) {}
230 Real initialGuess1(
const Array &acf, Real mean_r2,
233 Real A4 = acf[0] + mean_r2*mean_r2;
235 Real A = mean_r2*mean_r2/A4;
238 Real gammaLower = A <= 1./3. - tol_level ? std::sqrt((1 - 3*A)/(3 - 3*A)) + tol_level :
Real(tol_level);
239 Garch11Constraint constraints(gammaLower, 1.0 - tol_level);
241 Real gamma = gammaLower + (1 - gammaLower) * 0.5;
242 beta = std::min(gamma, std::max(gamma * (1 - A) - B, 0.0));
244 omega = mean_r2 * (1 - gamma);
247 gamma = std::max(gammaLower, -(1+4*B*B)/(4*B));
248 beta = std::min(gamma, std::max(gamma * (1 - A) - B, 0.0));
250 omega = mean_r2 * (1 - gamma);
253 gamma = std::max(gammaLower, -(1+B*B)/(2*B));
254 beta = std::min(gamma, std::max(gamma * (1 - A) - B, 0.0));
256 omega = mean_r2 * (1 - gamma);
258 Real D = (3*A-1)*(2*B*B+(1-A)*(2*A-1));
260 Real d = std::sqrt(D);
263 if (
b >= tol_level &&
b <= 1.0 - tol_level) {
264 g = (
b + B) / (1 - A);
266 if (g < gammaLower) {
268 if (
b >= tol_level &&
b <= 1.0 - tol_level) {
269 g = (
b + B) / (1 - A);
272 if (g >= gammaLower) {
274 beta = std::min(gamma, std::max(gamma * (1 - A) - B, 0.0));
276 omega = mean_r2 * (1 - gamma);
282 std::vector<std::size_t> idx;
283 std::size_t nCov = acf.size() - 1;
284 for (std::size_t i = 0; i <= nCov; ++i) {
285 if (i < 2 || (acf[i] > 0 && acf[i-1] > 0 && acf[i-1] > acf[i])) {
295 FitAcfConstraint c(gammaLower, 1.0 - tol_level);
296 NonLinearLeastSquare nnls(c);
297 nnls.setInitialValue(x);
298 FitAcfProblem pr(mean_r2, acf, idx);
299 x = nnls.perform(pr);
301 guess[0] = mean_r2 * (1 - x[0]);
302 guess[1] = x[0] - x[1];
304 if (constraints.test(guess)) {
309 }
catch (
const std::exception &) {
318 Real initialGuess2 (
const Array &acf, Real mean_r2,
321 Real A4 = acf[0] + mean_r2*mean_r2;
322 Real A = mean_r2*mean_r2/A4;
324 Real gammaLower = A <= 1./3. - tol_level ? std::sqrt((1 - 3*A)/(3 - 3*A)) + tol_level :
Real(tol_level);
325 Garch11Constraint constraints(gammaLower, 1.0 - tol_level);
330 std::vector<std::size_t> idx;
331 std::size_t nCov = acf.size() - 1;
332 for (std::size_t i = 0; i <= nCov; ++i) {
333 if (i < 2) idx.push_back(i);
334 if (i > 1 && acf[i] > 0 && acf[i-1] > 0 && acf[i-1] > acf[i]) {
335 gamma += acf[i]/acf[i-1];
342 if (gamma < gammaLower) gamma = gammaLower;
343 beta = std::min(gamma, std::max(gamma * (1 - A) - B, 0.0));
344 omega = mean_r2 * (1 - gamma);
351 FitAcfConstraint c(gammaLower, 1 - tol_level);
352 NonLinearLeastSquare nnls(c);
353 nnls.setInitialValue(x);
354 FitAcfProblem pr(mean_r2, acf, idx);
355 x = nnls.perform(pr);
357 guess[0] = mean_r2 * (1 - x[0]);
358 guess[1] = x[0] - x[1];
360 if (constraints.test(guess)) {
365 }
catch (
const std::exception &) {
377 auto cur = quoteSeries.
cbegin();
378 Real u = cur->second;
380 while (++cur != quoteSeries.
end()) {
381 sigma2 = omega +
alpha * u * u +
beta * sigma2;
382 retval[cur->first] = std::sqrt(sigma2);
385 sigma2 = omega +
alpha * u * u +
beta * sigma2;
388 retval[cur->first + (cur->first - (--prev)->first) ] = std::sqrt(sigma2);
393 ext::shared_ptr<Problem> Garch11::calibrate_r2(
394 Mode mode,
const std::vector<Volatility> &r2,
Real mean_r2,
396 EndCriteria endCriteria(10000, 500, tol_level, tol_level, tol_level);
398 return calibrate_r2(mode, r2, mean_r2, method, endCriteria,
402 ext::shared_ptr<Problem> Garch11::calibrate_r2(
403 Mode mode,
const std::vector<Volatility> &r2,
Real mean_r2,
411 "Data series is too short to fit GARCH model");
412 QL_REQUIRE (mean_r2 > 0,
"Data series is constant");
413 omega = mean_r2 * dataSize / (dataSize - 1);
416 Size maxLag = (
Size)std::sqrt(dataSize);
418 std::vector<Volatility> tmp(r2.size());
419 std::transform (r2.begin(), r2.end(), tmp.begin(),
420 [=](
Real x) ->
Real { return x - mean_r2; });
422 QL_REQUIRE (acf[0] > 0,
"Data series is constant");
424 Garch11CostFunction cost (r2);
427 Real gammaLower = 0.0;
430 if (mode != GammaGuess) {
431 gammaLower = initialGuess1(acf, mean_r2, opt1[1], opt1[2], opt1[0]);
432 fCost1 = cost.value(opt1);
437 if (mode != MomentMatchingGuess) {
438 gammaLower = initialGuess2(acf, mean_r2, opt2[1], opt2[2], opt2[0]);
439 fCost2 = cost.value(opt2);
442 Garch11Constraint constraints(gammaLower, 1.0 - tol_level);
444 ext::shared_ptr<Problem> ret;
445 if (mode != DoubleOptimization) {
447 ret = calibrate_r2(r2, method, constraints, endCriteria,
448 fCost1 <= fCost2 ? opt1 : opt2,
450 }
catch (
const std::exception &) {
451 if (fCost1 <= fCost2) {
462 ext::shared_ptr<Problem> ret1, ret2;
464 ret1 = calibrate_r2(r2, method, constraints, endCriteria,
469 if (constraints.test(opt1))
470 fCost1 = std::min(fCost1, cost.value(opt1));
471 }
catch (
const std::exception &) {
476 ret2 = calibrate_r2(r2, method, constraints, endCriteria,
481 if (constraints.test(opt2))
482 fCost2 = std::min(fCost2, cost.value(opt2));
483 }
catch (
const std::exception &) {
487 if (fCost1 <= fCost2) {
502 ext::shared_ptr<Problem> Garch11::calibrate_r2(
503 const std::vector<Volatility> &r2,
507 Garch11Constraint constraints(0.0, 1.0 - tol_level);
508 return calibrate_r2(r2, method, constraints, endCriteria,
512 ext::shared_ptr<Problem> Garch11::calibrate_r2(
513 const std::vector<Volatility> &r2,
518 std::vector<Volatility> tmp(r2.size());
519 std::transform(r2.begin(), r2.end(), tmp.begin(),
520 [=](
Real x) ->
Real { return x - mean_r2; });
521 return calibrate_r2(tmp, method, endCriteria, initGuess,
525 ext::shared_ptr<Problem> Garch11::calibrate_r2(
526 const std::vector<Volatility> &r2,
531 Garch11CostFunction cost(r2);
532 ext::shared_ptr<Problem> problem(
533 new Problem(cost, constraints, initGuess));
536 method.
minimize(*problem, endCriteria);
537 const Array &optimum = problem->currentValue();
544 ext::shared_ptr<Problem> Garch11::calibrate_r2(
545 const std::vector<Volatility> &r2,
551 std::vector<Volatility> tmp(r2.size());
552 std::transform(r2.begin(), r2.end(), tmp.begin(),
553 [=](
Real x) ->
Real { return x - mean_r2; });
554 return calibrate_r2(tmp, method, constraints, endCriteria,
autocovariance and convolution calculation
1-D array used in linear algebra.
const_iterator begin() const
Criteria to end optimization process:
Abstract class for constrained optimization method.
virtual EndCriteria::Type minimize(Problem &P, const EndCriteria &endCriteria)=0
minimize the optimization problem P
Constrained optimization problem.
Multi-dimensional simplex class.
Container for historical data.
const_iterator cbegin() const
const_iterator end() const
#define QL_REQUIRE(condition, message)
throw an error if the given pre-condition is not verified
ext::function< Real(Real)> b
std::vector< std::size_t > idx_
const std::vector< Volatility > & r2_
std::size_t Size
size of a container
Least square cost function.
void autocovariances(ForwardIterator begin, ForwardIterator end, OutputIterator out, std::size_t maxLag)
Unbiased auto-covariances.
Simplex optimization method.