<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/counterfactual-estimation-on-nonstationary-data-be-careful/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="https://insightr.wordpress.com/2017/07/31/counterfactual-estimation-on-nonstationary-data-be-careful/"> R – insightR</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p><strong>By Gabriel Vasconcelos</strong></p>
<p>In a recent paper that can be downloaded <a rel="nofollow" href="https://papers.ssrn.com/sol3/papers.cfm?abstract_id=2894065" target="_blank">here</a>, Carvalho, Masini and Medeiros show that estimating counterfactuals in a non-stationary framework (when I say non-stationary it means integrated) is a tricky task. It is intuitive that the models will not work properly in the absence of cointegration (spurious case), but what the authors show is that even with cointegration, the average treatment effect (ATE) converges to a non-standard distribution. As a result, standard tests on the ATE will identify treatment effects in several cases that there is no effect at all.</p>
<p><span></span></p>
<p>For those who are not familiar with counterfactual models: normally, these models have a treated unit (we want to know the effects of the treatment) and several untreated units that we call peers. There units may be cities, countries, companies, etc. Assuming that only one unit was treated and that there is no contamination, we can use information from the peers to project the behaviour of the treated unit as if it was not treated, which is the counterfactual. The models should also be able to identify if the treatment had no effects on the treated unit.</p>
<h2>Simulation tests</h2>
<p>Here I am going to generate cointegrated data with no treatment to show the behaviour of Synthetic Control and the ArCo (Other techniques will have the same features). The data is very simple and it is based on the mentioned article’s example. We are going to generate 100 observations of 10 time-series, 9 random walks and one that is the sum of these random walks plus an error, which is the treated unit. However, I will not include any treatment and I want the models to show that there was no treatment. The hypothesis to be tested is that at <img alt="t_0=51" title="t_0=51" src="https://s0.wp.com/latex.php?latex=t_0%3D51&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"> we had some treatment in the treated unit.</p>
<pre>library(ArCo)
library(glmnet)
library(Synth)

T=100 # Number of Observations
n=10 # Number of Time Series
set.seed(1) # Seed for replication

# Generate random walk peers
peers=matrix(0,T,n-1)
for(i in 2:T){
  peers[i,]=peers[i-1,]+rnorm(n-1)
}

# Generate the treated unit
y1=rowSums(peers)+rnorm(T)

# Put all the TS together
Y=cbind(y1,peers)

# Plot all TS. The black line is the &quot;treated unit&quot;
matplot(Y,type=&quot;l&quot;)
</pre>
<p><img alt="plot of chunk unnamed-chunk-16" src="https://i0.wp.com/i.imgur.com/FfxYo2B.png"></p>
<p>First, I am going to estimate the ArCo. The <code>$delta</code> show the average treatment effect with 95\% confidence interval. As you can see, the test indicated that the treatment was statistically significant (the zero is very far from the confidence interval). The plot also show a dislocated counterfactual (blue line).</p>
<pre># Estimate the ArCo in the non-stationary data
arco=fitArCo(list(Y),cv.glmnet,predict,treated.unity = 1,t0=51)
arco$delta
</pre>
<pre>##                  LB     delta        UB
## Variable1 -6.157713 -5.236705 -4.315697
</pre>
<pre>plot(arco)
</pre>
<p><img alt="plot of chunk unnamed-chunk-17" src="https://i2.wp.com/i.imgur.com/K210MOS.png"></p>
<p>Now I am going to estimate the Synthetic Control. The Synth package does not have an implementation of a formal test. However, we can have a good idea from the counterfactual plot, which also indicates that the treatment was effective (even though there was no treatment).</p>
<pre># Estimate the Synthetic Control in the nonstationary data

# Put the data in the Synth package notation
sY=as.vector(Y)
timeid=rep(1:T,ncol(Y))
unitid=sort(rep(1:ncol(Y),T))
synthdata=data.frame(time=timeid,unit=unitid,Y=sY)
synthdataprep&amp;lt;-
  dataprep(
    foo = synthdata,
    predictors = c(&quot;Y&quot;),
    predictors.op = &quot;mean&quot;,
    dependent = &quot;Y&quot;,
    unit.variable = &quot;unit&quot;,
    time.variable = &quot;time&quot;,
    treatment.identifier = 1,
    controls.identifier = c(2:10),
    time.predictors.prior = c(1:50),
    time.optimize.ssr = c(1:50),
    time.plot = 1:100
)
SC=synth(synthdataprep)
</pre>
<pre>path.plot(SC,synthdataprep)
abline(v=50,col=&quot;blue&quot;,lty=2)
</pre>
<p><img alt="plot of chunk unnamed-chunk-18" src="https://i0.wp.com/i.imgur.com/0Di7Cjj.png"></p>
<p>As you can see, both the ArCo and the Synthetic Control wrongly rejected the null hypothesis of no treatment. The best way to solve the problem is to take the first difference of all non-stationary time-series. However, there is an important thing to point out here. Imagine the random walk:</p>
<p><img alt="displaystyle y_t = y_{t-1} + varepsilon_t" title="displaystyle y_t = y_{t-1} + varepsilon_t" src="https://s0.wp.com/latex.php?latex=%5Cdisplaystyle+y_t+%3D+y_%7Bt-1%7D+%2B+%5Cvarepsilon_t&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"></p>
<p>Suppose we have a treatment in <img alt="t=t_0" title="t=t_0" src="https://s0.wp.com/latex.php?latex=t%3Dt_0&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"> that simply makes <img alt="y_{t_0} = c + y_{t_0-1}+varepsilon_{t_0}" title="y_{t_0} = c + y_{t_0-1}+varepsilon_{t_0}" src="https://s0.wp.com/latex.php?latex=y_%7Bt_0%7D+%3D+c+%2B+y_%7Bt_0-1%7D%2B%5Cvarepsilon_%7Bt_0%7D&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"> only for <img alt="t=t_0" title="t=t_0" src="https://s0.wp.com/latex.php?latex=t%3Dt_0&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex">. For any <img alt="tneq t_0" title="tneq t_0" src="https://s0.wp.com/latex.php?latex=t%5Cneq+t_0&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"> we will have <img alt="Delta y_t=varepsilon_t" title="Delta y_t=varepsilon_t" src="https://s0.wp.com/latex.php?latex=%5CDelta+y_t%3D%5Cvarepsilon_t&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"> and for <img alt="t=t_0" title="t=t_0" src="https://s0.wp.com/latex.php?latex=t%3Dt_0&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"> we will have <img alt="Delta y_{t_0}=c+varepsilon_{t_0}" title="Delta y_{t_0}=c+varepsilon_{t_0}" src="https://s0.wp.com/latex.php?latex=%5CDelta+y_%7Bt_0%7D%3Dc%2B%5Cvarepsilon_%7Bt_0%7D&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex">. Therefore, if we take the first difference of <img alt="y_t" title="y_t" src="https://s0.wp.com/latex.php?latex=y_t&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"> the treatment will have an impact only at <img alt="t=t_0" title="t=t_0" src="https://s0.wp.com/latex.php?latex=t%3Dt_0&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex">, which makes it impossible for any counterfactual model to identify an intervention. In this example, we would need an intervention that changed the drift of the Random Walk in a way that for all <img alt="tgeq t_0" title="tgeq t_0" src="https://s0.wp.com/latex.php?latex=t%5Cgeq+t_0&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"> we have:</p>
<p><img alt="displaystyle y_t=d+y_{t-1}+varepsilon_t" title="displaystyle y_t=d+y_{t-1}+varepsilon_t" src="https://s0.wp.com/latex.php?latex=%5Cdisplaystyle+y_t%3Dd%2By_%7Bt-1%7D%2B%5Cvarepsilon_t&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"></p>
<p>This will be better illustrated in an example. The plot below shows the difference between the two treatments. The red line is the case we can’t identify if we take the first difference. The Blue line changes the drift of the Random Walk, and therefore we can identify. <strong>It is impossible, with the tools I am using here, to identify if we are in the black line or in the red line</strong>.</p>
<pre>y=rep(0,T)
yt1=yt2=y
d=1
c=10
for(i in 2:T){
 e=rnorm(1)
 y[i]=y[i-1]+e
 if(i==51){
   yt1[i]=c+yt1[i-1]+e
 }else{
   yt1[i]=yt1[i-1]+e
 }
 if(i&amp;gt;=51){
   yt2[i]=d+yt2[i-1]+e
 }else{
   yt2[i]=yt2[i-1]+e
 }
}
plot(yt2,type=&quot;l&quot;,col=4,ylab=&quot;&quot;)
lines(yt1,col=2)
lines(y,col=1)
legend(&quot;topleft&quot;,legend=c(&quot;untreated&quot;,&quot;treated by a constant in t0&quot;,&quot;treated in the drift&quot;),col=c(1,2,4),bty=&quot;n&quot;,seg.len=0.5,cex=0.8,lty=1)
</pre>
<p><img alt="plot of chunk unnamed-chunk-19" src="https://i2.wp.com/i.imgur.com/KB25y8T.png"></p>
<p>Now let’s get back to the first example. The results clearly indicated an intervention where there was none. Now I will show what happens if we take the first difference.</p>
<pre># Take the first difference
dY=diff(Y,1)
# Estimate the ArCo
darco=fitArCo(list(dY),cv.glmnet,predict,treated.unity = 1,t0=50)
darco$delta
</pre>
<pre>##                   LB      delta        UB
## Variable1 -0.6892162 0.02737802 0.7439722
</pre>
<pre>plot(darco)
</pre>
<p><img alt="plot of chunk unnamed-chunk-20" src="https://i1.wp.com/i.imgur.com/eZjqJxr.png"></p>
<pre># Adjust the data to the Synth and estimate the model
dsY=as.vector(dY)
timeid=rep(2:T,ncol(Y))
unitid=sort(rep(1:ncol(Y),T-1))
dsynthdata=data.frame(time=timeid,unit=unitid,Y=dsY)

dsynthdataprep&amp;lt;-
  dataprep(
    foo = dsynthdata,
    predictors = c(&quot;Y&quot;),
    predictors.op = &quot;mean&quot;,
    dependent = &quot;Y&quot;,
    unit.variable = &quot;unit&quot;,
    time.variable = &quot;time&quot;,
    treatment.identifier = 1,
    controls.identifier = c(2:10),
    time.predictors.prior = c(2:50),
    time.optimize.ssr = c(2:50),
    time.plot = 2:100
)

dSC=synth(dsynthdataprep)
</pre>
<pre>##
## X1, X0, Z1, Z0 all come directly from dataprep object.
##
##
## ****************
##  optimization over w weights: computing synthtic control unit
##
##
##
## ****************
## ****************
## ****************
##
## MSPE (LOSS V): 8.176593
##
## solution.v:
##  1
##
## solution.w:
##  2.9e-09 1.4e-08 3.7e-09 0.9999999 7.4e-09 2.9e-09 1.16e-08 6.9e-09 4e-09
</pre>
<pre>path.plot(dSC,dsynthdataprep)
abline(v=50,col=&quot;blue&quot;,lty=2)
</pre>
<p><img alt="plot of chunk unnamed-chunk-21" src="https://i2.wp.com/i.imgur.com/aO8bVAG.png"></p>
<p>As you can see, the <code>$delta</code> in the ArCo showed that the treatment effect was statistically zero, which is true. The Synthetic Control plot also showed that there was no difference cause by the intervention even though the model adjustment was not very good.</p>
<p>Finally, note that we have two scenarios. If there is no intervention and we do not take the first difference we will probably wrongly reject the null hypothesis of no intervention. However, if we are in the “red line” case of the plot and we take the first difference, we will not reject the null when it may be false depending on the type of variable we are dealing with.</p>
<h2>References</h2>
<p><a rel="nofollow" href="https://papers.ssrn.com/sol3/papers.cfm?abstract_id=2894065" target="_blank">Carvalho, Carlos, Ricardo Masini, and Marcelo C. Medeiros. “The perils of counterfactual analysis with integrated processes.” (2016).</a></p>
<p> <a rel="nofollow" href="http://feeds.wordpress.com/1.0/gocomments/insightr.wordpress.com/489/" target="_blank"><img border="0" alt="" src="http://feeds.wordpress.com/1.0/comments/insightr.wordpress.com/489/"></a> <img border="0" alt="" src="https://i1.wp.com/pixel.wp.com/b.gif?resize=1%2C1&ssl=1" data-recalc-dims="1"></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://insightr.wordpress.com/2017/07/31/counterfactual-estimation-on-nonstationary-data-be-careful/"> R – insightR</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/R6ykWVvoLAo" height="1"></a></div>