How do I test if a given policy was successful? I have some data for medical R&amp;D and sales, my professor has asked me to check if a policy implemented in 2016 has made any effects on R&amp;D and R&amp;D/Sales data.
I tried to plot growth rates over years, but he says we don't need growth rates, and it has to be solved just by mean R&amp;D expenditures each year. (Mean over multiple companies). How do I go about significance testing of policy change?
Let us assume a balanced panel of companies $i=1,...,N$ observed in time periods $t=1,...,T$.
Let $s_{it}$ denote R&amp;D expenses for company $i$ time $t$ and let
$$\bar s_t :=\frac{1}{N} \sum_i s_{it},$$
denote the time $t$ average.
Consider now a simple regression model
$$(1)\ \ s_{it} = \beta_0 + \beta_1I[t\geq 2016] + e_{it}$$
where $I[t\geq 2016]$ is indicator for the time being after policy implementation. I assume the policy is intended to stimulate R&amp;D investments by companies and hence we assume $\beta_1&gt;0$ if the policy has effect.
Using the normal equations for POOLED OLS estimator it is easy to show that
$$\hat \beta_0 = \frac{1}{NT_1} \sum_i \sum_t s_{it}I[t&lt; 2016],$$
where $T_1$ is number of time periods before treatment. Hence $\hat \beta_1$ is roughly the average R&amp;D expense level BEFORE treatment/policy implementation.
If model (1) is true then $plim \ \ \hat \beta_0 = \beta_0 = \mathbb E[s_{it} \lvert t&lt;2016]$ which is the expected level of R&amp;D expenses for any company for any year before treatment.
You can also show that
$$\hat \beta_1 = \frac{1}{NT_2} \sum_i s_{it} I[t \geq 2016]-\frac{1}{NT_1} \sum_i \sum_t s_{it}I[t&lt; 2016]$$
If model (1) is true then $plim \ \ \hat \beta_1 = \beta_1 = \mathbb E[s_{it} \lvert t\geq 2016]- \mathbb E[s_{it} \lvert t&lt;2016]$ which is the expected level of R&amp;D expenses for any company for any year AFTER treatment.
You could plot the series $\{s_t\}_{t=1}^T$ and get - if model (1) is true - something like this

However something else may be happening over time and while you do not have any control group of companies not affected by the policy implementation you can still potentially do better than estimate model (1).
Consider the following plot

where in each time period going from $s_t$ to $s_{t+1}$ there is a small increase in R&amp;D let us call it $\lambda$. Then in the year of treatment the increase is slightly smaller.
The blue and read line are found using the estimators above and it is clear that because they do not take the time trend into consideration they overestimate the effect. The data used for the plot was generated after the model
$$s_{it} = \beta_0 + \lambda \cdot t + \beta_1 I[t\geq 2016] + e_{it},$$
which includes the time trend.
Here in this model the effect is immediate. But the effect could also happen in the years after the treatment instead of only in the year of treatment hence you could have a model like
$$s_{it} = \beta_0 + \lambda \cdot t + \beta_1 I[t\geq 2016] + \beta_2 I[
t\geq 2017]+ e_{it}.$$
The take away is that you should do as the professor says and calculate the time averages $s_t$ and plot them and look at them. Paying particular attention to what goes on around in the year treatment and in the years close to but still after treatment (based on the principle that effect comes after cause and takes some time be fully realized)
I R to simulate data and generate the plot by means of the following lines of code:
# simulate panel of RD
library(data.table)
set.seed(1)
N &lt;- 50 #Number of firms
T &lt;- 16     #Number of time periods
TT &lt;- 9 #policy implementation year (treatment year)

b_0 &lt;- 5
b_1 &lt;- 2

# Simulate model (1)
s &lt;- b_0 + b_1 * (1:T &gt;= TT) + rnorm(N*T)

# Make dataset 
id &lt;- rep(1:N,each=T)
time &lt;- rep(1:T,N)
I &lt;- as.numeric(time&gt;=TT)
mydata &lt;- data.table(id,time,s,I)

# Compute time averages
time_agg &lt;- mydata[,.(s_t=mean(s)),by=time]

#Plot time averages
plot(time_agg$s_t,ylim=c(0,10),pch=20,xlab=&quot;time&quot;,ylab=&quot;R&amp;D&quot;)
points(1:TT,rep(b_0,TT),type=&quot;l&quot;,lwd=1,col=&quot;blue&quot;)
points((TT):T,rep(b_0+b_1,T-TT+1),type=&quot;l&quot;,lwd=1,col=&quot;red&quot;)
savePlot(&quot;plot1.jpg&quot;,type=&quot;jpg&quot;)





# Simulate model (2)
s &lt;- b_0 + b_1 * (1:T &gt;= TT) + 0.5*1:T + rnorm(N*T)

# Make dataset 
mydata &lt;- data.table(id,time,s,I)

# Compute time averages
time_agg &lt;- mydata[,.(s_t=mean(s)),by=time]
model &lt;- lm(s~I,data=mydata)
hat_b_0 &lt;- coef(model)[1]
hat_b_1 &lt;- coef(model)[2]

#Plot time averages
plot(time_agg$s_t,ylim=c(5,16),pch=20,xlab=&quot;time&quot;,ylab=&quot;R&amp;D&quot;)
points(1:TT,rep(hat_b_0,TT),type=&quot;l&quot;,lwd=1,col=&quot;blue&quot;)
points((TT):T,rep(hat_b_0+hat_b_1,T-TT+1),type=&quot;l&quot;,lwd=1,col=&quot;red&quot;)
savePlot(&quot;plot2.jpg&quot;,type=&quot;jpg&quot;)