1/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3/*
4 Copyright (C) 2010 Liquidnet Holdings, Inc.
5
6 This file is part of QuantLib, a free-software/open-source library
7 for financial quantitative analysts and developers - http://quantlib.org/
8
9 QuantLib is free software: you can redistribute it and/or modify it
10 under the terms of the QuantLib license. You should have received a
11 copy of the license along with this program; if not, please email
12 <quantlib-dev@lists.sf.net>. The license is also available online at
13 <http://quantlib.org/license.shtml>.
14
15 This program is distributed in the hope that it will be useful, but WITHOUT
16 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
17 FOR A PARTICULAR PURPOSE. See the license for more details.
18*/
19
20/*! \file autocovariance.hpp
21 \brief autocovariance and convolution calculation
22*/
23
24#ifndef quantlib_auto_covariance_hpp
25#define quantlib_auto_covariance_hpp
26
27#include <ql/math/fastfouriertransform.hpp>
28#include <ql/math/array.hpp>
29#include <complex>
30#include <vector>
31#include <algorithm>
32#include <functional>
33
34namespace QuantLib {
35
36 namespace detail {
37
38 // Outputs double FT for a given input:
39 // input -> FFT -> norm -> FFT -> out
40 template <typename ForwardIterator>
41 std::vector<std::complex<Real> > double_ft(ForwardIterator begin,
42 ForwardIterator end) {
43 std::size_t nData = std::distance(begin, end);
44 std::size_t order = FastFourierTransform::min_order(inputSize: nData) + 1;
45 FastFourierTransform fft(order);
46 std::vector<std::complex<Real> > ft(fft.output_size());
47 fft.transform(begin, end, ft.begin());
48 Array tmp(ft.size(), 0.0);
49 std::complex<Real> z = std::complex<Real>();
50 for (Size i=0; i<ft.size(); ++i) {
51 tmp[i] = std::norm(z: ft[i]);
52 ft[i] = z;
53 }
54 fft.transform(inBegin: tmp.begin(), inEnd: tmp.end(), out: ft.begin());
55 return ft;
56 }
57
58
59 // Calculates and subtracts mean from the input data; returns mean
60 template <typename InputIterator, typename OutputIterator>
61 Real remove_mean(InputIterator begin, InputIterator end,
62 OutputIterator out) {
63 Real mean(0.0);
64 std::size_t n = 1;
65 for (InputIterator it = begin; it != end; ++it, ++n)
66 mean = (mean*Real(n-1) + *it)/n;
67 std::transform(begin, end, out, [=](Real x) -> Real { return x - mean; });
68 return mean;
69 }
70
71 }
72
73
74 //! Convolutions of the input sequence.
75 /*! Calculates x[0]*x[n]+x[1]*x[n+1]+x[2]*x[n+2]+...
76 for n = 0,1,...,maxLag via FFT.
77
78 \pre The size of the output sequence must be maxLag + 1
79 */
80 template <typename ForwardIterator, typename OutputIterator>
81 void convolutions(ForwardIterator begin, ForwardIterator end,
82 OutputIterator out, std::size_t maxLag) {
83 using namespace detail;
84 std::size_t nData = std::distance(begin, end);
85 QL_REQUIRE(maxLag < nData, "maxLag must be less than data size");
86 const std::vector<std::complex<Real> >& ft = double_ft(begin, end);
87 Real w = 1.0 / (Real)ft.size();
88 for (std::size_t k = 0; k <= maxLag; ++k)
89 *out++ = ft[k].real() * w;
90 }
91
92 //! Unbiased auto-covariances
93 /*! Results are calculated via FFT.
94
95 \pre Input data are supposed to be centered (i.e., zero mean).
96 \pre The size of the output sequence must be maxLag + 1
97 */
98 template <typename ForwardIterator, typename OutputIterator>
99 void autocovariances(ForwardIterator begin, ForwardIterator end,
100 OutputIterator out, std::size_t maxLag) {
101 using namespace detail;
102 std::size_t nData = std::distance(begin, end);
103 QL_REQUIRE(maxLag < nData,
104 "number of covariances must be less than data size");
105 const std::vector<std::complex<Real> >& ft = double_ft(begin, end);
106 Real w1 = 1.0 / (Real)ft.size(), w2 = (Real)nData;
107 for (std::size_t k = 0; k <= maxLag; ++k, w2 -= 1.0) {
108 *out++ = ft[k].real() * w1 / w2;
109 }
110 }
111
112 //! Unbiased auto-covariances
113 /*! Results are calculated via FFT.
114
115 This overload accepts non-centered data, removes the mean and
116 returns it as a result. The centered sequence is written back
117 into the input sequence if the reuse parameter is true.
118
119 \pre The size of the output sequence must be maxLag + 1
120 */
121 template <typename ForwardIterator, typename OutputIterator>
122 Real autocovariances(ForwardIterator begin, ForwardIterator end,
123 OutputIterator out,
124 std::size_t maxLag, bool reuse) {
125 using namespace detail;
126 Real mean = 0.0;
127 if (reuse) {
128 mean = remove_mean(begin, end, begin);
129 autocovariances(begin, end, out, maxLag);
130 } else {
131 Array tmp(std::distance(begin, end));
132 mean = remove_mean(begin, end, tmp.begin());
133 autocovariances(tmp.begin(), tmp.end(), out, maxLag);
134 }
135 return mean;
136 }
137
138
139 //! Unbiased auto-correlations.
140 /*! Results are calculated via FFT.
141 The first element of the output is the unbiased sample variance.
142
143 \pre Input data are supposed to be centered (i.e., zero mean).
144 \pre The size of the output sequence must be maxLag + 1
145 */
146 template <typename ForwardIterator, typename OutputIterator>
147 void autocorrelations(ForwardIterator begin, ForwardIterator end,
148 OutputIterator out, std::size_t maxLag) {
149 using namespace detail;
150 std::size_t nData = std::distance(begin, end);
151 QL_REQUIRE(maxLag < nData,
152 "number of correlations must be less than data size");
153 const std::vector<std::complex<Real> >& ft = double_ft(begin, end);
154 Real w1 = 1.0 / (Real)ft.size(), w2 = (Real)nData;
155 Real variance = ft[0].real() * w1 / w2;
156 *out++ = variance * w2 / (w2-1.0);
157 w2 -= 1.0;
158 for (std::size_t k = 1; k <= maxLag; ++k, w2 -= 1.0)
159 *out++ = ft[k].real() * w1 / (variance * w2);
160 }
161
162 //! Unbiased auto-correlations.
163 /*! Results are calculated via FFT.
164 The first element of the output is the unbiased sample variance.
165
166 This overload accepts non-centered data, removes the mean and
167 returns it as a result. The centered sequence is written back
168 into the input sequence if the reuse parameter is true.
169
170 \pre The size of the output sequence must be maxLag + 1
171 */
172 template <typename ForwardIterator, typename OutputIterator>
173 Real autocorrelations(ForwardIterator begin, ForwardIterator end,
174 OutputIterator out,
175 std::size_t maxLag, bool reuse) {
176 using namespace detail;
177 Real mean = 0.0;
178 if (reuse) {
179 mean = remove_mean(begin, end, begin);
180 autocorrelations(begin, end, out, maxLag);
181 } else {
182 Array tmp(std::distance(begin, end));
183 mean = remove_mean(begin, end, tmp.begin());
184 autocorrelations(tmp.begin(), tmp.end(), out, maxLag);
185 }
186 return mean;
187 }
188
189}
190
191#endif
192

source code of quantlib/ql/math/autocovariance.hpp