1 // (C) Copyright Nick Thompson 2019.
2 // Use, modification and distribution are subject to the
3 // Boost Software License, Version 1.0. (See accompanying file
4 // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
5
6 #ifndef BOOST_MATH_STATISTICS_LJUNG_BOX_HPP
7 #define BOOST_MATH_STATISTICS_LJUNG_BOX_HPP
8
9 #include <cmath>
10 #include <iterator>
11 #include <utility>
12 #include <boost/math/distributions/chi_squared.hpp>
13 #include <boost/math/statistics/univariate_statistics.hpp>
14
15 namespace boost::math::statistics {
16
17 template<class RandomAccessIterator>
ljung_box(RandomAccessIterator begin,RandomAccessIterator end,int64_t lags=-1,int64_t fit_dof=0)18 auto ljung_box(RandomAccessIterator begin, RandomAccessIterator end, int64_t lags = -1, int64_t fit_dof = 0) {
19 using Real = typename std::iterator_traits<RandomAccessIterator>::value_type;
20 int64_t n = std::distance(begin, end);
21 if (lags >= n) {
22 throw std::domain_error("Number of lags must be < number of elements in array.");
23 }
24
25 if (lags == -1) {
26 // This is the same default as Mathematica; it seems sensible enough . . .
27 lags = static_cast<int64_t>(std::ceil(std::log(Real(n))));
28 }
29
30 if (lags <= 0) {
31 throw std::domain_error("Must have at least one lag.");
32 }
33
34 auto mu = boost::math::statistics::mean(begin, end);
35
36 std::vector<Real> r(lags + 1, Real(0));
37 for (size_t i = 0; i < r.size(); ++i) {
38 for (auto it = begin + i; it != end; ++it) {
39 Real ak = *(it) - mu;
40 Real akml = *(it-i) - mu;
41 r[i] += ak*akml;
42 }
43 }
44
45 Real Q = 0;
46
47 for (size_t k = 1; k < r.size(); ++k) {
48 Q += r[k]*r[k]/(r[0]*r[0]*(n-k));
49 }
50 Q *= n*(n+2);
51
52 typedef boost::math::policies::policy<
53 boost::math::policies::promote_float<false>,
54 boost::math::policies::promote_double<false> >
55 no_promote_policy;
56
57 auto chi = boost::math::chi_squared_distribution<Real, no_promote_policy>(Real(lags - fit_dof));
58
59 Real pvalue = 1 - boost::math::cdf(chi, Q);
60 return std::make_pair(Q, pvalue);
61 }
62
63
64 template<class RandomAccessContainer>
ljung_box(RandomAccessContainer const & v,int64_t lags=-1,int64_t fit_dof=0)65 auto ljung_box(RandomAccessContainer const & v, int64_t lags = -1, int64_t fit_dof = 0) {
66 return ljung_box(v.begin(), v.end(), lags, fit_dof);
67 }
68
69 }
70 #endif
71