<p><div><div><div><div data-show-faces="false" data-height="21" data-width="100" data-layout="button_count" data-send="true" data-href="https://www.r-bloggers.com/time-series-and-mcht/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="https://ntguardian.wordpress.com/2018/11/12/time-series-and-mcht/"> R – Curtis Miller's Personal Website</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<h2>Introduction</h2>
<p>Over the past few weeks I’ve published articles about my new package, <strong>MCHT</strong>, starting with an <a rel="nofollow" href="https://ntguardian.wordpress.com/2018/10/08/announcing-mcht-package-bootstrap-monte-carlo-hypothesis-testing/" target="_blank">introduction</a>, a <a rel="nofollow" href="https://wp.me/p27liZ-TG" target="_blank">further technical discussion</a>, demonstrating <a rel="nofollow" href="https://wp.me/p27liZ-TQ" target="_blank">maximized Monte Carlo (MMC) hypothesis testing</a>, <a rel="nofollow" href="https://wp.me/p27liZ-TU" target="_blank">bootstrap hypothesis testing</a>, and last week I showed <a rel="nofollow" href="https://wp.me/p27liZ-U0" target="_blank">how to handle multi-sample and multivariate data</a>. This is the final article where I explain the capabilities of the package. I show how <strong>MCHT</strong> can handle time series data.</p>
<p><span></span></p>
<p>I should mention that I’m not focused on the merits of the procedures I use as examples in these posts, and that’s going to be the case here. It’s possible (perhaps even likely) that there’s a better way to decide between the hypotheses than what I show here. In these articles, I’m more interested in showing what <em>can</em> be done rather than what <em>should</em> be done. In particular, I like simple examples that many can understand, even if they may not be the best tool for the task at hand.</p>
<p>So far I don’t think this has been a serious issue; that is, I don’t think the procedures I’ve shown so far could be considered controversial (I think the most controversial would be the permutation test example). But the example I want to use here could be argued with; I personally would not use it. That said, I’m still willing to demonstrate it because it doesn’t take much to understand what’s going on and it does demonstrate how time series data can be handled.</p>
<h2>Context</h2>
<p>Suppose we want to perform a test for the location of the mean, and thus decide between the hypotheses</p>
<p><img alt="H_0: mu = mu_0" title="H_0: mu = mu_0" src="https://s0.wp.com/latex.php?latex=H_0%3A+%5Cmu+%3D+%5Cmu_0&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"></p>
<p><img alt="H_A: mu neq mu_0" title="H_A: mu neq mu_0" src="https://s0.wp.com/latex.php?latex=H_A%3A+%5Cmu+%5Cneq+%5Cmu_0&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"></p>
<p>There is the usual <img alt="t" title="t" src="https://s0.wp.com/latex.php?latex=t&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex">-statistic, which is <img alt="frac{bar{X} - mu_0}{frac{S}{sqrt{n}}}" title="frac{bar{X} - mu_0}{frac{S}{sqrt{n}}}" src="https://s0.wp.com/latex.php?latex=%5Cfrac%7B%5Cbar%7BX%7D+-+%5Cmu_0%7D%7B%5Cfrac%7BS%7D%7B%5Csqrt%7Bn%7D%7D%7D&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex">, and as mentioned before the statistic assumes that the data came from a Normal distribution. That’s not all the test assumes, though. It also assumes that the data is independent and identically distributed.</p>
<p>In cross-sectional contexts this is fine, but it’s not okay when the data could depend on time and thus is not independent and identically distributed. Suppose instead that our data was generated according to a <a rel="nofollow" href="https://en.wikipedia.org/wiki/Autoregressive%E2%80%93moving-average_model" target="_blank">first-order autoregressive process (AR(1))</a>, described below:</p>
<p><img alt="X_t = mu + epsilon_t" title="X_t = mu + epsilon_t" src="https://s0.wp.com/latex.php?latex=X_t+%3D+%5Cmu+%2B+%5Cepsilon_t&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"></p>
<p><img alt="epsilon_t = rho epsilon_{t - 1} + sigma w_t" title="epsilon_t = rho epsilon_{t - 1} + sigma w_t" src="https://s0.wp.com/latex.php?latex=%5Cepsilon_t+%3D+%5Crho+%5Cepsilon_%7Bt+-+1%7D+%2B+%5Csigma+w_t&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"></p>
<p>In this context, assume <img alt="w_t sim N(0, 1)" title="w_t sim N(0, 1)" src="https://s0.wp.com/latex.php?latex=w_t+%5Csim+N%280%2C+1%29&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> and is independent and identically distribution. It’s no longer given that the conventional <img alt="t" title="t" src="https://s0.wp.com/latex.php?latex=t&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex">-test will work as marketed since the data is no longer independent or identically distributed. Additionally, we have two nuisance parameters, <img alt="rho" title="rho" src="https://s0.wp.com/latex.php?latex=%5Crho&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> and <img alt="sigma" title="sigma" src="https://s0.wp.com/latex.php?latex=%5Csigma&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex">, that need to be accounted for.</p>
<h2>Monte Carlo Testing with Time Series</h2>
<p>We will view <img alt="rho" title="rho" src="https://s0.wp.com/latex.php?latex=%5Crho&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> and <img alt="sigma" title="sigma" src="https://s0.wp.com/latex.php?latex=%5Csigma&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> as nuisance parameters and use MMC testing to handle them. That leaves the question of how to simulate an AR(1) process. With <strong>MCHT</strong>, if you can simulate a process, you can test with it.</p>
<p>The <a rel="nofollow" href="https://en.wikipedia.org/wiki/Time_series" target="_blank">time series model above has a stationary solution</a> when <img alt="|rho| &lt; 1" title="|rho| &lt; 1" src="https://s0.wp.com/latex.php?latex=%7C%5Crho%7C+%3C+1&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> and when <img alt="t" title="t" src="https://s0.wp.com/latex.php?latex=t&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> ranges between <img alt="-infty" title="-infty" src="https://s0.wp.com/latex.php?latex=-%5Cinfty&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> and <img alt="infty" title="infty" src="https://s0.wp.com/latex.php?latex=%5Cinfty&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex">. It's not possible to simulate a series of infinite length but one can get close by simulating a series that is very long. In particular, one can simulate, say, 500 terms of the series starting at a fixed number, then the actual number of terms of the series wanted, then throw away the first 500 terms. This is known as burn-in and it's very common practice in time series simulation.</p>
<p>Fortunately <code>MCHTest()</code> allows for burn-in. Suppose that the sample size of the actual dataset is <img alt="n" title="n" src="https://s0.wp.com/latex.php?latex=n&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> and we've decided that we want a burn-in period of <img alt="B" title="B" src="https://s0.wp.com/latex.php?latex=B&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex">. Then we can do the following:</p>
<ol>
<li>Generate <img alt="n + B" title="n + B" src="https://s0.wp.com/latex.php?latex=n+%2B+B&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> random numbers to represent <img alt="w_t" title="w_t" src="https://s0.wp.com/latex.php?latex=w_t&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> (except possibly for the scaling factor, as we're treating that as a nuisance parameter).</li>
<li>Apply the recursive formula described above to the <img alt="w_t" title="w_t" src="https://s0.wp.com/latex.php?latex=w_t&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> series after scaling the series by <img alt="sigma" title="sigma" src="https://s0.wp.com/latex.php?latex=%5Csigma&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> and using a chosen <img alt="rho" title="rho" src="https://s0.wp.com/latex.php?latex=%5Crho&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex">, and add <img alt="mu" title="mu" src="https://s0.wp.com/latex.php?latex=%5Cmu&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> to it.</li>
<li>Keep only the last <img alt="n" title="n" src="https://s0.wp.com/latex.php?latex=n&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> terms of the series; throw away the rest. This is your simulated dataset.</li>
<li>After having obtained the simulated dataset, proceed with the Monte Carlo test as usual.</li>
</ol>
<p>With MMC, the unscaled series <img alt="w_t" title="w_t" src="https://s0.wp.com/latex.php?latex=w_t&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> is fixed after we generate it and we use optimization to adversarially choose <img alt="rho" title="rho" src="https://s0.wp.com/latex.php?latex=%5Crho&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> and <img alt="sigma" title="sigma" src="https://s0.wp.com/latex.php?latex=%5Csigma&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> so that we maximize the <img alt="p" title="p" src="https://s0.wp.com/latex.php?latex=p&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex">-value of the test.</p>
<p>When using <code>MCHTest()</code>, the <code>rand_gen</code> function does not need to produce a dataset of the same length as the original dataset; this allows for burning it. However, if you're going to do this, then the <code>stat_gen</code> function needs to know what the sample size of the dataset is, but all you need to do is give the <code>stat_gen</code> function the parameter <code>n</code>; this will be given the sample size of the original dataset. And of course the <code>test_stat</code> function won't care whether the data came from a time series or not.</p>
<p>Putting this all together, we create the following test.</p>
<pre>
library(MCHT)
library(doParallel)

registerDoParallel(detectCores())

ts &lt;- function(x, mu = 0) {
  sqrt(length(x)) * (mean(x) - mu)/sd(x)
}

rg &lt;- function(n) {
  rnorm(n + 500)  # Extra terms for a burn-in period
}

sg &lt;- function(x, n, mu = 0, rho = 0, sigma = 1) {
  x &lt;- sigma * x
  if (abs(rho) >= 1) {stop(&quot;Bad rho given!&quot;)}
  eps &lt;- filter(x, rho, &quot;recursive&quot;)  # Apply the recursion
  eps &lt;- eps[-(1:500)]  # Throw away first 500 observations; they're burn-in
  dat &lt;- eps + mu

  test_stat(dat, mu = mu)  # Will be localizing
}

mc.ar1.t.test &lt;- MCHTest(ts, sg, rg, N = 1000, seed = 123, test_params = &quot;mu&quot;,
                         nuisance_params = c(&quot;rho&quot;, &quot;sigma&quot;),
                         optim_control = list(lower = c(&quot;rho&quot; = -0.999,
                                                        &quot;sigma&quot; = 0),
                                              upper = c(&quot;rho&quot; = 0.999,
                                                        &quot;sigma&quot; = 100),
                                              control = list(&quot;max.time&quot; = 10)),
                         threshold_pval = 0.2, localize_functions = TRUE,
                         lock_alternative = FALSE)

dat &lt;- c(-1.02, -1.13, 0.53, 0.21, 1.76, 1.79, 1.42, -0.31, -0.28, -0.44)

mc.ar1.t.test(dat, mu = 0, alternative = &quot;two.sided&quot;)
</pre>
<pre>
## Warning in mc.ar1.t.test(dat, mu = 0, alternative = &quot;two.sided&quot;): Computed
## p-value is greater than threshold value (0.2); the optimization algorithm
## may have terminated early
</pre>
<pre>
## 
## 	Monte Carlo Test
## 
## data:  dat
## S = 0.73415, p-value = 0.264
## alternative hypothesis: true mu is not equal to 0
</pre>
<pre>
mc.ar1.t.test(dat, mu = 3, alternative = &quot;two.sided&quot;)
</pre>
<pre>
## Warning in mc.ar1.t.test(dat, mu = 3, alternative = &quot;two.sided&quot;): Computed
## p-value is greater than threshold value (0.2); the optimization algorithm
## may have terminated early
</pre>
<pre>
## 
## 	Monte Carlo Test
## 
## data:  dat
## S = -7.9712, p-value = 0.504
## alternative hypothesis: true mu is not equal to 3
</pre>
<pre>
t.test(dat, mu = 3, alternative = &quot;two.sided&quot;)  # For reference
</pre>
<pre>
## 
## 	One Sample t-test
## 
## data:  dat
## t = -7.9712, df = 9, p-value = 2.278e-05
## alternative hypothesis: true mean is not equal to 3
## 95 percent confidence interval:
##  -0.5265753  1.0325753
## sample estimates:
## mean of x 
##     0.253
</pre>
<h2>Conclusion</h2>
<p>I have now covered what I consider the essential technical functionality of <strong>MCHT</strong>. All of the functionality I described in these posts is functionality that I want this package to have. Thus I personally am quite happy this package exists, which is good; I'm the package's primary audience, after all. All I can hope is that others find the package useful too.</p>
<p>I wrote this article more than a month before it was published, so perhaps I have made an update that isn't being accounted for here, but as of this version (0.1.0), I'd call the package in a beta stage of stability; it's usable, but features could be added or removed and there could be unknown bugs.</p>
<p>The following is a list of possible areas of expansion. This list exists mostly because I think it needs to exist; it gives me something to aim for before making a 1.0 release. That said, they could be useful features.</p>
<p>*A function for making diagnostic-type plots for tests, such as a function creating a plot for the rejection probability function (RPF) as described in [1].<br>
*A function that accepts a <code>MCHTest</code>-class object and returns a function that, rather than returning a <code>htest</code>-class object, returns a function that will give the test statistic, simulated test statistics, and a <img alt="p" title="p" src="https://s0.wp.com/latex.php?latex=p&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex">-value, in a list; could be useful for diagnostic work.<br>
*Real-world datasets that can be used for examples.<br>
*Functions with a simpler interface than <code>MCHTest</code>, perhaps with more restrictions on inputs.<br>
*Pre-made <code>MCHTest</code> objects perhaps implementing common Monte Carlo or Bootstrap tests.</p>
<p>I also welcome community requests and collaboration. If you want a feature, consider issuing a pull request on GitHub.</p>
<p>Do you want more documentation? More examples? More background? Let me know! I'd be willing to write more on this subject. Perhaps if I amass enough content I could write a book documenting <strong>MCHT</strong> and Monte Carlo/bootstrap testing.</p>
<p>These blog posts together extend beyond 10,000 words, so I'm thinking I have enough material to submit an article to, say, <a rel="nofollow" href="https://www.jstatsoft.org/" target="_blank"><em>J. Stat. Soft.</em></a> or the <a rel="nofollow" href="https://journal.r-project.org/" target="_blank"><em>R Journal</em></a> and thus get my first publication where I'm the sole author. But this is something I'm still considering; I'm an insecure person at heart.</p>
<p>Next week I will still be using this package in a blog post, but I won't be writing about how to use it anymore; instead, I'll be using it to revisit <a rel="nofollow" href="https://ntguardian.wordpress.com/2018/04/19/a-recession-before-2020-is-likely-on-the-distribution-of-time-between-recessions/" target="_blank">a proposition I made many months ago</a>. (It was because of that article I created this package.) Stay tuned, and thanks for reading!</p>
<h2>References</h2>
<ol>
<li>R. Davidson and J. G. MacKinnon, <em>The size distortion of bootstrap test</em>, Econometric Theory, vol. 15 (1999) pp. 361-376</li>
</ol>
<hr>
<p><a rel="nofollow" href="https://www.packtpub.com/" target="_blank">Packt Publishing</a> published a book for me entitled <a rel="nofollow" href="https://ntguardian.wordpress.com/books/hands-on-data-analysis-with-numpy-and-pandas/" target="_blank"><em>Hands-On Data Analysis with NumPy and Pandas</em></a>, a book based on my video course <a rel="nofollow" href="https://ntguardian.wordpress.com/video-courses/unpacking-numpy-and-pandas/" target="_blank"><em>Unpacking NumPy and Pandas</em></a>. This book covers the basics of setting up a Python environment for data analysis with Anaconda, using Jupyter notebooks, and using NumPy and pandas. If you are starting out using Python for data analysis or know someone who is, please consider <a rel="nofollow" href="https://www.packtpub.com/big-data-and-business-intelligence/hands-data-analysis-numpy-and-pandas" target="_blank">buying my book</a> or at least spreading the word about it. You can buy the book directly or purchase a <a rel="nofollow" href="https://mapt.packtpub.com/?gclid=Cj0KCQjwuuHdBRCvARIsAELQRQGrXDoR4ox216bUlSD6mTQYibN24CA5bYF7QikiQQ6j3cGm63HRXsEaAjJjEALw_wcB" target="_blank">subscription to Mapt</a> and read it there.</p>
<p>If you like my blog and would like to support it, spread the word (if not get a copy yourself)!</p>
	<p>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="https://ntguardian.wordpress.com/2018/11/12/time-series-and-mcht/"> R – Curtis Miller's Personal Website</a></strong>.</div>
<hr>
<a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="https://feedburner.google.com/fb/a/mailverify?uri=RBloggers">daily e-mail updates</a></strong> about <a rel="nofollow" href="https://www.r-project.org/" title="The R Project for Statistical Computing">R</a> news and <a rel="nofollow" href="https://www.r-bloggers.com/search/tutorial" title="R tutorials">tutorials</a> on topics such as: <a rel="nofollow" href="https://www.r-bloggers.com/search/data%20science" title="Data science">Data science</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Big%20Data" title="Big Data">Big Data, <a rel="nofollow" href="https://www.r-users.com/" title="R jobs">R jobs</a>, visualization (<a rel="nofollow" href="https://www.r-bloggers.com/search/ggplot2" title="ggplot and ggplot2 tutorials">ggplot2</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/boxplot" title="Boxplots using lattice and ggplot2 tutorials">Boxplots</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/map" title="Maps and gis">maps</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/animation" title="Animation in R">animation</a>), programming (<a rel="nofollow" href="https://www.r-bloggers.com/search/RStudio" title="RStudio IDE for R">RStudio</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/sweave" title="Sweave and literate programming">Sweave</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX" title="LaTeX in R">LaTeX</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/SQL" title="SQL and databases">SQL</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/eclipse" title="Eclipse IDE for R">Eclipse</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/git" title="git and github, Version Control System">git</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/hadoop" title="Large data in R using Hadoop">hadoop</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping" title="Web Scraping of google, facebook, yahoo, twitter and more using R">Web Scraping</a>) statistics (<a rel="nofollow" href="https://www.r-bloggers.com/search/regression" title="Regressions and ANOVA analysis tutorials">regression</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/PCA" title="principal component analysis tutorial">PCA</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/time+series" title="Time series">time series</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/trading" title="finance trading">trading</a>) and more...
<img width="1" alt="" src="http://feeds.feedburner.com/~r/RBloggers/~4/Kpnr_KyPH4Q" height="1"></a></div>