boost/math/statistics/ljung_box.hpp

   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>
  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>
  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