71 lines · plain
1// (C) Copyright Nick Thompson 2019.2// Use, modification and distribution are subject to the3// Boost Software License, Version 1.0. (See accompanying file4// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)5 6#ifndef BOOST_MATH_STATISTICS_LJUNG_BOX_HPP7#define BOOST_MATH_STATISTICS_LJUNG_BOX_HPP8 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 15namespace boost::math::statistics {16 17template<class RandomAccessIterator>18auto 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 64template<class RandomAccessContainer>65auto 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#endif71