<p><div><div><div><div data-width="100" data-show-faces="false" data-send="true" data-layout="button_count" data-height="21" data-href="https://www.r-bloggers.com/can-we-predict-stock-prices-with-prophet/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="http://msperlin.github.io/2017-03-05-Prophet-and_stock-market/"> R and Finance</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p> – </p>
<p>Facebook recently released a API package allowing access to its<br>
forecasting model called<br>
<a target="_blank" rel="nofollow" href="http://blog.revolutionanalytics.com/2017/02/facebook-prophet.html">prophet</a>.<br>
According to the underling post:</p>
<div>
<pre><code>It's not your traditional ARIMA-style time series model. It's closer in spirit to a  Bayesian-influenced generalized additive model, a regression of smooth terms. The model is resistant   to the effects of outliers, and supports data collected over an irregular time scale (ingliding presence of missing data) without the need for interpolation. The underlying calculation engine is Stan; the R and Python packages simply provide a convenient interface.  
</code></pre>
</div>
<p>After reading it, I got really curious about the predictive performance<br>
of this method for stock prices. That is, <strong>can we predict stock price<br>
movements based on prophet?</strong> In this post I will investigate this<br>
research question using a database of prices for the SP500 components.</p>
<p>Before describing the code and results, it is noteworthy to point out<br>
that forecasting stock returns is really hard! There is a significant<br>
body of literature trying to forecast prices and to prove (or not) that<br>
financial markets are efficient in pricing publicly available<br>
information, including historical prices. This is the so called<br>
efficient market hypothesis. I have studied it, tried to trade for<br>
myself for a while when I was a Msc student, advised several graduate<br>
students on it, and the results are mostly the same: it is very<br>
difficult to find a trade signal that works well and is sustainable in<br>
real life.</p>
<p>This means that most of the variation in prices is due to random factors<br>
that cannot be anticipated. The explanation is simple, prices move<br>
according to investor’s expectation from available information. Every<br>
time that new (random) information, true or not, reaches the market,<br>
investor’s update their beliefs and trade accordingly. So, unless, new<br>
information or market expectation have a particular pattern, price<br>
changes will be mostly random.</p>
<p>Even with a body of evidence against our research, it is still<br>
interesting to see how we could apply <code>prophet</code> in a trading setup.</p>
<h2>The data</h2>
<p>First, let’s download stock prices for all components of the SP500 index<br>
in the previous two years. This is a lengthy download. If you don’t have<br>
the patience or time, I saved the the result in RData file available<br>
<a target="_blank" rel="nofollow" href="http://msperlin.github.io/data/SP500.RData">here</a>.</p>
<div>
<pre><code>library(BatchGetSymbols)

my.stocks &lt;- GetSP500Stocks()$ticker

first.date &lt;- as.Date('2015-01-01')
last.date &lt;- Sys.Date()
df.stocks &lt;- BatchGetSymbols(my.stocks, 
                             first.date = first.date, 
                             last.date = last.date)[[2]]

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
</code></pre>
</div>
<p>In order to make the code run faster, I will select only 25 random<br>
stocks from the sample:</p>
<div>
<pre><code>set.seed(10)
n.stocks &lt;- 25
my.stocks &lt;- sample(unique(df.stocks$ticker), n.stocks)

df.stocks &lt;- df.stocks[df.stocks$ticker %in% my.stocks, ]
</code></pre>
</div>
<p>Now, let’s understand how prophet works. I was happy to see that the<br>
interface is quite simple, you offer a time series with input <code>y</code> and a<br>
date vector with <code>ds</code>. If no further custom option is set, you are good<br>
to go. My only complain with <code>prophet</code> is that that the function outputs<br>
lots of messages. They really should add a <code>quiet</code> option, so that the<br>
user doesn’t have to use <code>capture.output</code> to silent it. Have a look in<br>
the next example with a dummy series:</p>
<div>
<pre><code>library(prophet)

## Loading required package: Rcpp

df.est &lt;- data.frame(y = rnorm(100), ds = Sys.Date() + 1:100)

m &lt;- prophet(df = df.est)

## STAN OPTIMIZATION COMMAND (LBFGS)
## init = user
## save_iterations = 1
## init_alpha = 0.001
## tol_obj = 1e-012
## tol_grad = 1e-008
## tol_param = 1e-008
## tol_rel_obj = 10000
## tol_rel_grad = 1e+007
## history_size = 5
## seed = 492642109
## initial log joint probability = -8.50656
## Optimization terminated normally: 
##   Convergence detected: absolute parameter change was below tolerance
</code></pre>
</div>
<p>The next step is to think about how to structure a function for our<br>
research problem. Our study has two steps, first we will set a<br>
<em>training</em> (in-sample) period, estimate the model and make forecasts.<br>
After that, we use the <em>out-of-sample</em> data to test the accuracy of the<br>
model.</p>
<p>The whole procedure of estimating and forecasting will be encapsulated<br>
in a single R function. This is not the best way of doing it but, for<br>
our simple example, it will suffice. My function will take as input a<br>
dataframe and the number of out-of-sample forecasts. Based on the<br>
adjusted closing prices, we calculate returns and feed<br>
<code>1:(nrow(df)-nfor)</code> rows for the estimation. The last <code>nfor</code> rows are<br>
used for testing the accuracy of the model. For example, if I have a<br>
vector with 1000 returns and <code>nfor=5</code>, I use observations from <code>1:995</code><br>
for estimating the model and <code>996:1000</code> for testing the forecasts. The<br>
function returns a dataframe with the predictions for each horizon, its<br>
error, among other things. Here’s the function definition:</p>
<div>
<pre><code>est.model.and.forecast &lt;- function(df.in, nfor=5){
  # Estimated model using prophet and forecast it
  #
  # Args:
  #   df.in - A dataframe with columns price.adjusted and ref.date
  #   nfor - Number of out-of-sample forecasts
  #
  # Returns:
  #   A dataframe with forecasts and errors for each horizon.
  
  require(prophet)
  require(dplyr)
  
  my.ticker &lt;- df.in$ticker[1]
  
  #cat('\nProcessing ', my.ticker)
  
  n.row &lt;- nrow(df.in)
  df.in$ret &lt;- with(df.in, c(NA,price.adjusted[2:n.row]/price.adjusted[1:(n.row - 1)] - 1))
  
  df.in &lt;- select(df.in, ref.date, ret)
  names(df.in) &lt;- c('ds', 'y')
  
  idx &lt;- nrow(df.in) - nfor
  
  df.est &lt;- df.in[1:idx, ]
  df.for &lt;- df.in[(idx + 1):nrow(df.in), ]
  
  capture.output(
    m &lt;- prophet(df = df.est)
  )
  
  # forecast 50 days ahead (it also includes non trading days)
  df.pred &lt;- predict(m,
                     make_future_dataframe(m,
                                           periods = nfor + 50))
  
  
  df.for &lt;- merge(df.for, df.pred, by = 'ds')
  df.for &lt;- select(df.for, ds, y, yhat)
  
  # forecast statistics
  df.for$eps &lt;- with(df.for,y - yhat)
  df.for$abs.eps &lt;- with(df.for,abs(y - yhat))
  df.for$perc.eps &lt;- with(df.for,(y - yhat)/y)
  df.for$nfor &lt;- 1:nrow(df.for)
  df.for$ticker &lt;- my.ticker
  
  return(df.for)
  
}
</code></pre>
</div>
<p>Let’s try it out using the <code>by</code> function to apply it for each stock in<br>
our sample. All results are later combined in a single dataframe with<br>
function <code>do.call</code>.</p>
<div>
<pre><code>out.l &lt;- by(data = df.stocks,
            INDICES = df.stocks$ticker, 
            FUN = est.model.and.forecast, nfor = 5)

my.result &lt;- do.call(rbind, out.l)
</code></pre>
</div>
<p>Lets have a look in the resulting dataframe:</p>
<div>
<pre><code>head(my.result)

##                ds            y          yhat           eps      abs.eps
## AN.1   2017-02-27  0.004697779  1.073144e-03  0.0036246355 0.0036246355
## AN.2   2017-02-28 -0.024442020 -7.927652e-04 -0.0236492544 0.0236492544
## AN.3   2017-03-01  0.015468387  2.439663e-03  0.0130287240 0.0130287240
## AN.4   2017-03-02 -0.003647329  6.270727e-05 -0.0037100361 0.0037100361
## AN.5   2017-03-03 -0.013350539 -6.233167e-04 -0.0127272219 0.0127272219
## ANTM.1 2017-02-27  0.006592211  6.073838e-03  0.0005183736 0.0005183736
##          perc.eps nfor ticker
## AN.1   0.77156364    1     AN
## AN.2   0.96756548    2     AN
## AN.3   0.84228070    3     AN
## AN.4   1.01719266    4     AN
## AN.5   0.95331149    5     AN
## ANTM.1 0.07863425    1   ANTM
</code></pre>
</div>
<p>In this object you’ll find the forecasts (yhat), the actual values (y),<br>
the absolute and normalized error (abs.eps, perc.eps).</p>
<p>For ou first analysis, let’s have a look on the effect of the<br>
forecasting horizon over the absolute error distribution.</p>
<div>
<pre><code>library(ggplot2)

p &lt;- ggplot(my.result, aes(x=factor(nfor), 
                           y=abs.eps))
p &lt;- p + geom_boxplot()

print(p)
</code></pre>
</div>
<p><img alt="" src="https://i2.wp.com/msperlin.github.io/img/2017-03-05-Prophet-and_stock-market_files/figure-markdown_strict/unnamed-chunk-8-1.png?w=456" data-recalc-dims="1"></p>
<p>We do find some positive dependency. As the horizon increases, the<br>
forecasting algorithm makes more mistakes. Surprisingly, this pattern is<br>
not found for <code>nfor=5</code>. It might be interesting to add more data and<br>
check if this effect is robust.</p>
<h2>Encopassing test</h2>
<p>A simple and powerful test for verifying the accuracy of a prediction<br>
algorithm is the encompassing test. The idea is to estimate the<br>
following linear model with the real returns (<em>R</em><sub><em>t</em></sub>) and its<br>
predictions ($\hat{R} _t$).</p>
<p><img alt="" src="https://i1.wp.com/msperlin.github.io/2017-03-05-Prophet-and_stock-market/EncompassingTest.png?w=456" data-recalc-dims="1"></p>
<p>If the model provides good forecasts, we can expect that <em>α</em> is equal to<br>
zero (no bias) and <em>β</em> is equal to 1. If both conditions are true, we<br>
have that $R_t = \hat{R} _t + \epsilon _t$$, meaning that our<br>
forecasting model provides an unbiased estimator of the predicted<br>
variable. In a formal research, we could use a Wald test to verify this<br>
hypothesis jointly.</p>
<p>First, lets find the result of the encompassing test for all forecasts.</p>
<div>
<pre><code>lm.model &lt;- lm(formula = y ~yhat, data = my.result)
summary(lm.model)

## 
## Call:
## lm(formula = y ~ yhat, data = my.result)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.041310 -0.005200 -0.000846  0.006056  0.035647 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)   
## (Intercept) -0.002240   0.001289  -1.738  0.08479 . 
## yhat         0.822045   0.254414   3.231  0.00158 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01133 on 123 degrees of freedom
## Multiple R-squared:  0.07824,    Adjusted R-squared:  0.07074 
## F-statistic: 10.44 on 1 and 123 DF,  p-value: 0.001582
</code></pre>
</div>
<p>As you can see, it didn’t work very well. While the constant seems Ok,<br>
the value of 0.8220446 is not even close to 1. But, it could be the case<br>
that the different horizon have different results. A longer horizon,<br>
with bad forecasts, will be affecting short horizons with good<br>
forecasts. Lets use <code>dplyr</code> to separate our model according to <code>nfor</code>.</p>
<div>
<pre><code>models &lt;- my.result %>%
  group_by(nfor) %>%
  do(ols.model = lm(data = ., formula = y ~ yhat ))
</code></pre>
</div>
<p>We report the results with <code>broom::tidy</code>.</p>
<div>
<pre><code>library(broom)
knitr::kable(models %>% tidy(ols.model))
</code></pre>
</div>
<table>
<thead>
<tr>
<th align="right">nfor</th>
<th align="left">term</th>
<th align="right">estimate</th>
<th align="right">std.error</th>
<th align="right">statistic</th>
<th align="right">p.value</th>
</tr>
</thead>
<tbody>
<tr>
<td align="right">1</td>
<td align="left">(Intercept)</td>
<td align="right">0.0024381</td>
<td align="right">0.0017693</td>
<td align="right">1.3779940</td>
<td align="right">0.1814651</td>
</tr>
<tr>
<td align="right">1</td>
<td align="left">yhat</td>
<td align="right">-0.3374193</td>
<td align="right">0.4488191</td>
<td align="right">-0.7517936</td>
<td align="right">0.4597999</td>
</tr>
<tr>
<td align="right">2</td>
<td align="left">(Intercept)</td>
<td align="right">-0.0071114</td>
<td align="right">0.0022989</td>
<td align="right">-3.0933634</td>
<td align="right">0.0051276</td>
</tr>
<tr>
<td align="right">2</td>
<td align="left">yhat</td>
<td align="right">1.2326600</td>
<td align="right">0.4826284</td>
<td align="right">2.5540562</td>
<td align="right">0.0177381</td>
</tr>
<tr>
<td align="right">3</td>
<td align="left">(Intercept)</td>
<td align="right">0.0074791</td>
<td align="right">0.0026961</td>
<td align="right">2.7740540</td>
<td align="right">0.0107937</td>
</tr>
<tr>
<td align="right">3</td>
<td align="left">yhat</td>
<td align="right">1.1998427</td>
<td align="right">0.4618100</td>
<td align="right">2.5981305</td>
<td align="right">0.0160770</td>
</tr>
<tr>
<td align="right">4</td>
<td align="left">(Intercept)</td>
<td align="right">-0.0079042</td>
<td align="right">0.0017917</td>
<td align="right">-4.4114616</td>
<td align="right">0.0002019</td>
</tr>
<tr>
<td align="right">4</td>
<td align="left">yhat</td>
<td align="right">0.3345002</td>
<td align="right">0.3219130</td>
<td align="right">1.0391012</td>
<td align="right">0.3095587</td>
</tr>
<tr>
<td align="right">5</td>
<td align="left">(Intercept)</td>
<td align="right">-0.0041302</td>
<td align="right">0.0030427</td>
<td align="right">-1.3574059</td>
<td align="right">0.1878234</td>
</tr>
<tr>
<td align="right">5</td>
<td align="left">yhat</td>
<td align="right">0.7356029</td>
<td align="right">0.6071194</td>
<td align="right">1.2116281</td>
<td align="right">0.2379564</td>
</tr>
</tbody>
</table>
<p>Well, it again didn’t work as expected. We do not find that shorter<br>
horizons have better results in the encompassing test. In fact, we even<br>
get some negative betas! This means that, for some horizons, it might be<br>
better to take the opposite suggestion of the forecast!</p>
<h2>Trading based on forecasts</h2>
<p>In a practical trading applications, it might not be of interest to<br>
forecast actual returns. If you are trading according to these<br>
forecasts, you are probably more worried about the direction of the<br>
forecasts and not its nominal error. A model can have bad nominal<br>
forecasts, but be good in predicting the sign of the next price<br>
movement. If this is the case, you can still make money even though your<br>
model fails in the encompassing test.</p>
<p>Let’s try it out with a simple trading strategy for all different<br>
horizons:</p>
<ul>
<li>buy in end of day <em>t</em> if forecast in <em>t+1</em> is positive and sell at<br>
the end of <em>t+1</em></li>
<li>short-sell in the end of day <em>t</em> when forecast for <em>t+1</em> is negative<br>
and buy it back in the end of <em>t+1</em></li>
</ul>
<p>The total profit will be given by:</p>
<div>
<pre><code>my.profit &lt;- sum(with(my.result, (yhat>0)*y + (yhat&lt;0)*-y))
print(my.profit)

## [1] 0.07947572
</code></pre>
</div>
<p>Not bad! Doesn’t look like much, but remember that we have a few trading<br>
days and this return might be due to a sistematic effect in the market.<br>
Let’s see how this result compares to random trading signals.</p>
<div>
<pre><code>n.sim &lt;- 10000

monkey.ret &lt;- numeric(length = n.sim)
for (i in seq(n.sim)) {
  rnd.vec &lt;- rnorm(length(my.result$y))
  
  monkey.ret[i] &lt;- sum( (rnd.vec>0)*my.result$y + (rnd.vec&lt;0)*-my.result$y )
  
} 

temp.df &lt;- data.frame(monkey.ret, my.profit)
p &lt;- ggplot(temp.df, aes(monkey.ret)) 
p &lt;- p + geom_histogram()
p &lt;- p + geom_vline(aes(xintercept =  my.profit),size=2)
p &lt;- p + labs(x='Returns from random trading signals')
print(p)

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
</code></pre>
</div>
<p><img alt="" src="https://i0.wp.com/msperlin.github.io/img/2017-03-05-Prophet-and_stock-market_files/figure-markdown_strict/unnamed-chunk-13-1.png?w=456" data-recalc-dims="1"></p>
<p>The previous histogram shows the total return from randomnly generated<br>
signals in 10^{4} simulations. The vertical line is the result from<br>
using <code>prophet</code>. As you can see, it is a bit higher than the average of<br>
the distribution. The total return from <code>prophet</code> is lower than the<br>
return of the naive strategy in 27.5 percent of the simulations. This is<br>
not a bad result. But, notice that we didnt add trading or liquidity<br>
costs to the analysis, which will make the total returns worse.</p>
<p>The main results of this simple study are clear: <strong><code>prophet</code> is bad at<br>
point forecasts for returns but does quite better in directional<br>
predictions</strong>. It might be interesting to test it further, with more<br>
data, adding trading costs, other forecasting setups, and see if the<br>
results hold.</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="http://msperlin.github.io/2017-03-05-Prophet-and_stock-market/"> R and Finance</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 title="The R Project for Statistical Computing" rel="nofollow" href="https://www.r-project.org/">R</a> news and <a title="R tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/tutorial">tutorials</a> on topics such as: <a title="Data science" rel="nofollow" href="https://www.r-bloggers.com/search/data%20science">Data science</a>, <a title="Big Data" rel="nofollow" href="https://www.r-bloggers.com/search/Big%20Data">Big Data, <a title="R jobs" rel="nofollow" href="https://www.r-users.com/">R jobs</a>, visualization (<a title="ggplot and ggplot2 tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/ggplot2">ggplot2</a>, <a title="Boxplots using lattice and ggplot2 tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/boxplot">Boxplots</a>, <a title="Maps and gis" rel="nofollow" href="https://www.r-bloggers.com/search/map">maps</a>, <a title="Animation in R" rel="nofollow" href="https://www.r-bloggers.com/search/animation">animation</a>), programming (<a title="RStudio IDE for R" rel="nofollow" href="https://www.r-bloggers.com/search/RStudio">RStudio</a>, <a title="Sweave and literate programming" rel="nofollow" href="https://www.r-bloggers.com/search/sweave">Sweave</a>, <a title="LaTeX in R" rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX">LaTeX</a>, <a title="SQL and databases" rel="nofollow" href="https://www.r-bloggers.com/search/SQL">SQL</a>, <a title="Eclipse IDE for R" rel="nofollow" href="https://www.r-bloggers.com/search/eclipse">Eclipse</a>, <a title="git and github, Version Control System" rel="nofollow" href="https://www.r-bloggers.com/search/git">git</a>, <a title="Large data in R using Hadoop" rel="nofollow" href="https://www.r-bloggers.com/search/hadoop">hadoop</a>, <a title="Web Scraping of google, facebook, yahoo, twitter and more using R" rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping">Web Scraping</a>) statistics (<a title="Regressions and ANOVA analysis tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/regression">regression</a>, <a title="principal component analysis tutorial" rel="nofollow" href="https://www.r-bloggers.com/search/PCA">PCA</a>, <a title="Time series" rel="nofollow" href="https://www.r-bloggers.com/search/time+series">time series</a>, <a title="finance trading" rel="nofollow" href="https://www.r-bloggers.com/search/trading">trading</a>) and more...
<img height="1" alt="" width="1" src="http://feeds.feedburner.com/~r/RBloggers/~4/3EMXKPrUv9U"></a></div>