diff --git a/include/boost/math/statistics/ljung_box.hpp b/include/boost/math/statistics/ljung_box.hpp
new file mode 100644
index 000000000..597fac2f6
--- /dev/null
+++ b/include/boost/math/statistics/ljung_box.hpp
@@ -0,0 +1,68 @@
+//  (C) Copyright Nick Thompson 2019.
+//  Use, modification and distribution are subject to the
+//  Boost Software License, Version 1.0. (See accompanying file
+//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+#ifndef BOOST_MATH_STATISTICS_LJUNG_BOX_HPP
+#define BOOST_MATH_STATISTICS_LJUNG_BOX_HPP
+
+#include <cmath>
+#include <iterator>
+#include <utility>
+#include <boost/math/distributions/chi_squared.hpp>
+#include <boost/math/statistics/univariate_statistics.hpp>
+
+namespace boost::math::statistics {
+
+template<class RandomAccessIterator>
+auto ljung_box(RandomAccessIterator begin, RandomAccessIterator end, int64_t lags = -1, int64_t fit_dof = 0) {
+    using Real = typename std::iterator_traits<RandomAccessIterator>::value_type;
+    int64_t n = std::distance(begin, end);
+    if (lags >= n) {
+      throw std::domain_error("Number of lags must be < number of elements in array.");
+    }
+
+    if (lags == -1) {
+      // This is the same default as Mathematica; it seems sensible enough . . .
+      lags = static_cast<int64_t>(std::ceil(std::log(Real(n))));
+      std::cout << "Number of lags = " << lags << "\n";
+    }
+
+    if (lags <= 0) {
+      throw std::domain_error("Must have at least one lag.");
+    }
+
+
+    std::vector<Real> r(lags + 1, Real(0));
+    for (size_t i = 0; i < r.size(); ++i) {
+      for (auto it = begin + i; it != end; ++it) {
+        Real ak = *(it);
+        Real akml = *(it-i);
+        r[i] += ak*akml;
+      }
+    }
+
+    Real Q = 0;
+    for (size_t i = 1; i < r.size(); ++i) {
+      r[i] /= r[0];
+      std::cout << "r[" << i << "] = " << r[i] << "\n";
+    }
+
+
+    for (size_t k = 1; k < r.size(); ++k) {
+      Q += r[k]*r[k]/((n-k));
+    }
+    std::cout << "Q/(n*(n+2)) = " << Q << "\n";
+    Q *= n*(n+2);
+    Real pvalue = -1;
+    return std::make_pair(Q, pvalue);
+}
+
+
+template<class RandomAccessContainer>
+auto ljung_box(RandomAccessContainer const & v, int64_t lags = -1, int64_t fit_dof = 0) {
+    return ljung_box(v.begin(), v.end(), lags, fit_dof);
+}
+
+}
+#endif
diff --git a/test/ljung_box_test.cpp b/test/ljung_box_test.cpp
new file mode 100644
index 000000000..e0843977e
--- /dev/null
+++ b/test/ljung_box_test.cpp
@@ -0,0 +1,44 @@
+/*
+ * Copyright Nick Thompson, 2019
+ * Use, modification and distribution are subject to the
+ * Boost Software License, Version 1.0. (See accompanying file
+ * LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+ */
+
+#include "math_unit_test.hpp"
+#include <numeric>
+#include <utility>
+#include <random>
+#include <boost/math/statistics/ljung_box.hpp>
+
+using boost::math::statistics::ljung_box;
+
+template<class Real>
+void test_trivial()
+{
+  std::vector<Real> v{1,2};
+
+  double expected_statistic = Real(32)/Real(25);
+
+  auto [computed_statistic, computed_pvalue] = ljung_box(v, 1);
+  CHECK_ULP_CLOSE(expected_statistic, computed_statistic, 10);
+
+}
+
+
+void test_agreement_with_mathematica()
+{
+  std::vector<double> v{0.7739928761039216,-0.4468259278452086,0.98287381303903,-0.3943029116201079,0.6569015496559457};
+  double expected_statistic = 10.2076093223439;
+
+  auto [computed_statistic, computed_pvalue] = ljung_box(v);
+  CHECK_ULP_CLOSE(expected_statistic, computed_statistic, 10);
+
+}
+
+int main()
+{
+    test_trivial<double>();
+    //test_agreement_with_mathematica();
+    return boost::math::test::report_errors();
+}