# Phase matching condition for enhancement of phase sensitivity in quantum
metrology

Jing Liu, Xiaoxing Jing, Xiaoguang Wang<EMAIL_ADDRESS>Zhejiang
Institute of Modern Physics, Department of Physics, Zhejiang University,
Hangzhou 310027, China

###### Abstract

We find a phase matching condition for enhancement of sensitivity in a Mach-
Zehnder interferometer illuminated by an arbitrary state in one input port and
an odd(even) state in the other port. Under this condition, the Fisher
information becomes maximal with respect to the relative phase of two modes
and the phase sensitivity is enhanced. For the case with photon losses, we
further find that the phase matching condition keeps unchanged with a coherent
state and a coherent superposition state as the input states.

###### pacs:

03.67.-a, 42.50.St, 42.50.Dv

_Introduction_.–With the development of quantum information theory 1 ; 1-1 ;
1-1.1 ; 1-2 ; 1-3  and quantum technology 2 ; 2-1 ; 2-2 , quantum metrology is
becoming more and more practical nowadays. Improving the precision of a
parameter is always the basic theme in quantum metrology. The fundamental
procedure of it is the parametrization process. There are three mainly methods
to perform this procedure: (i) unitary parametrization; (ii) channel
parametrization; (iii) accelerating parametrization. Unitary parametrization
is the most useful and well studied method because it is widely applied in the
phase estimation, the main task in quantum metrology. The precision
measurement of gravity, temperature, week magnetic strength and many other
parameters can be classified in the category of phase estimation. Channel
parametrization is also widely studied in the channel estimation Fujiwara ;
channel_est  for many years. With the help of the quantum technology, it may
be possible to proceed the parametrization in a accelerating way relative , in
which the relativistic effects cannot be neglected any more. In this paper, we
mainly focus on the category of phase estimation and study how to effectively
enhance the phase sensitivity.

In 1981, Caves Caves  found out that, for phase estimation, taking a high
intensity coherent state and a low intensity squeezed vacuum state as the
input states of a Mach-Zehnder interferometer, the precision can beat the
shot-noise limit(quantum standard limit), i.e., $1/\sqrt{N}$, where $N$ is the
total photon number in both modes. Since then, many protocols have been
proposed to fulfill the similar job, such as NOON state, entangled coherent
state Joo , two-mode squeezed state Anisimov , number squeezed states
Pezze2013  and so on. With high intensity, some of these states can even
theoretically achieve or surpass the Heisenberg limit, i.e., $1/N$. In the
pioneer work of Caves Caves , to enhance the precision, the phases of the two
input states need to satisfy a relation. This can be considered as a kind of
phase matching condition (PMC). Thus, it is reasonable to study if there is a
more general PMC for more general states to to enhance the phase sensitivity.
This is the major motivation of this paper. To depict the precision of a
parameter $\theta$, quantum Fisher information(QFI) is an available useful
concept because it describes the lower bound on the variance of the estimator
$\hat{\theta}$ due to the Cramér-Rao theorem: $\mathrm{Var}(\hat{\theta})\geq
1/(\nu F)$, Helstrom ; Holevo  where $\mathrm{Var}(\cdot)$ is the variance,
$\nu$ is the number of repeated experiments and $F$ is the QFI.

Figure 1: (Color online) Sketch of the Mach-Zehnder interferometer. The
transformation of the beam splitters are described by $\exp(\pm i\pi
J_{x}/2)$, and that of the phase shift is described by $\exp(i\theta J_{z})$.
The input states in port A and B are an arbitrary state and an even(odd)
state, respectively.

In this paper, we discuss a general scenario of a Mach-Zehnder interferometer.
In this interferometer, one of the input port is an _arbitrary_ state and the
other one is an _even_(odd) state. This scenario covers many important cases
including the famous protocol of Caves Caves , and a recent one proposed by
Pezzé and Smerzi Pezze2013 . We give an analytic expression of the QFI and
identify a PMC to optimize the parameter precision. Under this condition, the
QFI is only determined the the average photon numbers of the two modes and the
corresponding expectation values of the square of annihilation operators. We
then give two examples of our scenario. Further, for the case that photon
losses occur in both arms with the same transmission coefficients, and the
input state is a product of a coherent state and a coherent superposition
state, the analytic expression of QFI is provided. Based on this expression,
we prove that the PMC keeps unchanged for any transmission coefficients.

_Mach-Zehnder interferometer_.–Mach-Zehnder(MZ) interferometer is a well known
optical device in quantum metrology which is constructed with two beam
splitters and one or two phase shifts. The interferometer we consider here is
constructed with two 50:50 beam splitters and two phase shifts in both arms,
as shown in Fig. 1. The well used 50:50 beam splitter can be described by
Campos ; Sanders ; Yurke  $B_{x}=\exp(-i\frac{\pi}{2}J_{x})$, where
$J_{x}=\frac{1}{2}(a^{\dagger}b+b^{\dagger}a)$ is one of the operators in the
Schwinger representation of bosons. The others are
$J_{y}=\frac{1}{2i}(a^{\dagger}b-b^{\dagger}a)$ and
$J_{z}=\frac{1}{2}\left(a^{\dagger}a-b^{\dagger}b\right)$. Here $a$, $b$ are
annihilation operators for ports A and B, respectively. Operators $J_{x}$,
$J_{y}$, $J_{z}$ satisfy the commutation:
$\left[J_{i},J_{j}\right]=i\epsilon_{ijk}J_{k}$, with $\epsilon_{ijk}$ the so-
called Levi-Civita symbol. Usually, a phase shift can be expressed by a
unitary transformation $\exp[i\theta\hat{N}_{\mathrm{A(B)}}]$, where
$\hat{N}_{\mathrm{A}}=a^{\dagger}a$, $\hat{N}_{\mathrm{B}}=b^{\dagger}b$ and
$\theta$ is an unknown parameter. In our scenario, denoting $\theta$ as the
relative phase between the two arms, one can describe the transformation of
the total phase shift by the operator $P_{z}=\exp(i\theta J_{z})$. With above
devices we can construct a well used MZ interferometer, the transformation of
which can be written as $U_{\mathrm{mz}}=B_{x}P_{z}B_{x}^{\dagger}$. Through
some algebra, one can find that this transformation can be simplified as

$\displaystyle U_{\mathrm{MZ}}$ $\displaystyle=$
$\displaystyle\exp\left(-i\theta J_{y}\right),$ (1)

which is a rotation along $y$ direction.

_Phase matching condition for QFI_.–Quantum Fisher information(QFI) is a
central concept in quantum metrology, and it is defined as Helstrom ; Holevo
$F:=\mathrm{Tr}\left(\rho L^{2}\right)$, where $L$ is the so-called symmetric
logarithmic derivative determined by
$\partial_{\theta}\rho_{\theta}=\left(\rho_{\theta}L+L\rho_{\theta}\right)/2.$
We consider a separable input state
$\rho_{\mathrm{in}}=\rho_{\mathrm{A}}\otimes\rho_{\mathrm{B}}$. Here
$\rho_{\mathrm{A}}$ is an arbitrary state and
$\rho_{\mathrm{B}}=|\varphi\rangle\langle\varphi|$ is an even(odd) state.
Utilizing the spectral decomposition
$\rho_{\mathrm{A}}=\sum_{j=1}^{M}p_{j}|\psi_{j}\rangle\langle\psi_{j}|$, the
quantum Fisher information can be written as Liu ; Helstrom ; Holevo

$F=\sum_{j=1}^{M}4p_{j}\langle\phi_{j}|J_{y}^{2}|\phi_{j}\rangle-\sum_{j,j^{\prime}=1}^{M}\frac{8p_{j}p_{j^{\prime}}}{p_{j}+p_{j^{\prime}}}|\langle\phi_{j}|J_{y}|\phi_{j^{\prime}}\rangle|^{2},$
(2)

where $M$ is the dimension of the support of $\rho_{\mathrm{A}}$ and
$|\phi_{j}\rangle=|\psi_{j}\rangle\otimes|\varphi\rangle$ is an eigenstate of
$\rho_{\mathrm{in}}$ with eigenvalue $p_{j}$.

As the input state in port B is an even(odd) state in our scenario, it
satisfies $\langle\varphi|b|\varphi\rangle=0$. Then after some algebra, the
QFI can be expressed by

$\displaystyle F$ $\displaystyle=$ $\displaystyle
2\bar{n}_{\mathrm{A}}\bar{n}_{\mathrm{B}}+\bar{n}_{\mathrm{A}}+\bar{n}_{\mathrm{B}}-2\mathrm{Re}\left(\langle
a^{\dagger 2}\rangle\langle b^{2}\rangle\right),$ (3)

where $\bar{n}_{\mathrm{A}}=\mathrm{Tr}(\rho_{\mathrm{A}}a^{\dagger}a)$ and
$\bar{n}_{\mathrm{B}}=\langle\varphi|b^{\dagger}b|\varphi\rangle$ are average
photon numbers for mode $a$ and $b$, respectively. To optimize this Fisher
information, the input states need to satisfy the following PMC:

$|\mathrm{Arg}(\langle a^{2}\rangle)-\mathrm{Arg}(\langle b^{2}\rangle)|=\pi.$
(4)

Under this condition, the QFI becomes

$F_{\mathrm{m}}=2\bar{n}_{\mathrm{A}}\bar{n}_{\mathrm{B}}+\bar{n}_{\mathrm{A}}+\bar{n}_{\mathrm{B}}+2|\langle
a^{2}\rangle||\langle b^{2}\rangle|.$ (5)

One key feature of the above result is that the QFI is only determined by the
mean photon numbers and the expectations of $a^{2}$ and $b^{2}$. We emphasize
that one input state is arbitrary.

For the case that $|\langle a^{2}\rangle||\langle b^{2}\rangle|=0$, the phases
of the input states can be chosen arbitrarily. For example, if
$\rho_{\mathrm{B}}$ is a Fock state, i.e., $\rho_{\mathrm{B}}=|n\rangle\langle
n|$, then the QFI reduces to
$F_{\mathrm{m}}=2\bar{n}_{\mathrm{A}}\bar{n}_{\mathrm{B}}+\bar{n}_{\mathrm{A}}+\bar{n}_{\mathrm{B}}$,
which is the result in Ref. Pezze2013  and has been discussed in detail. In
the following, we give two examples of the PMC.

Example 1.–In this example, we choose $\rho_{\mathrm{A}}$ be a coherent state
$|\beta\rangle\langle\beta|$ and $\rho_{\mathrm{B}}$ a coherent superposition
state $|{\alpha}\rangle_{+}\langle\alpha|$, where
$|{\alpha}\rangle_{+}=N_{\alpha}(|\alpha\rangle+|-\alpha\rangle)$ with
$|\alpha\rangle$ also a coherent state and
$N_{\alpha}^{2}=1/(2+2e^{-2|\alpha|^{2}})$. It is easy to find that
${}_{+}\langle{\alpha}|b|{\alpha}\rangle_{+}=0$. Here, we denote
$\alpha=|\alpha|\exp(i\Phi_{\alpha})$, $\beta=|\beta|\exp(i\Phi_{\beta})$,
then from Eq. (4), the phase matching condition can be specifically written as

$|\Phi_{\alpha}-\Phi_{\beta}|=\frac{\pi}{2}.$ (6)

Under this condition, the QFI reduces to

$F_{\mathrm{m}}=2\bar{n}_{\mathrm{A}}\bar{n}_{\mathrm{B}}+\bar{n}_{\mathrm{A}}+\bar{n}_{\mathrm{B}}+2\bar{n}_{\mathrm{A}}|\alpha|^{2},$
(7)

where $\bar{n}_{\mathrm{A}}=|\beta|^{2}$ and
$\bar{n}_{\mathrm{B}}=|\alpha|^{2}\tanh|\alpha|^{2}$. As $\tanh|\alpha|^{2}$
is a monotonic function and very close to $1$ for $|\alpha|\geq 2$, then for
most value of $|\alpha|$, $\bar{n}_{\mathrm{B}}$ is equal to $|\alpha|^{2}$,
and the QFI reduces to
$F_{\mathrm{m}}=4\bar{n}_{\mathrm{A}}\bar{n}_{\mathrm{B}}+\bar{n}_{\mathrm{A}}+\bar{n}_{\mathrm{B}}$.
It is not difficult to obtain that $F_{\mathrm{m}}\leq N^{2}+N$, where
$N=\bar{n}_{\mathrm{A}}+\bar{n}_{\mathrm{B}}$ is the total photon number. The
equality can be achieved when $\bar{n}_{\mathrm{A}}=\bar{n}_{\mathrm{B}}$. For
the case that total photon number $N$ is fixed, this bound is the optimal
value of the QFI. Considering the PMC, the optimal value can be achieved when
$\beta=\pm i\alpha$. From the optimal value of $F_{\mathrm{m}}=N^{2}+N\geq
N^{2}$, one can find that, with high intensity, the pair of the coherent state
$|\pm i\alpha\rangle\langle\pm i\alpha|$ and the coherent superposition state
$|{\alpha}\rangle_{+}$ can surpass the Heisenberg limit. Ref. Joo  considered
the same input states in our scenario but with only one phase shift in one
arm. In their case, utilizing the similar analysis above, it is easy to find
that the input states $|\alpha\rangle$ and $|{\alpha}\rangle_{+}$ is the
optimal choice and can reach the maximum value of QFI.

Figure 2: (Color online) The variation of $F_{\mathrm{m}}/N^{2}$ with the
average particle number $\bar{n}_{\mathrm{A}}$ and $\bar{n}_{\mathrm{B}}$. The
input states here are a coherent state $|\beta\rangle\langle\beta|$ in port A
and a squeezed vacuum state $|\xi\rangle\langle\xi|$ in port B.

Example 2.–Another well known even state is the squeezed vacuum state, which
is defined as Scully  $|\xi\rangle=S(\xi)|0\rangle$, where the squeezing
operator reads $S(\xi)=\exp(\frac{1}{2}\xi^{*}b^{2}-\frac{1}{2}\xi
b^{\dagger^{2}})$ with $\xi=|\xi|\exp(i\Phi_{\xi})$. For convenience, we still
choose the input state in port A as a coherent state
$|\beta\rangle\langle\beta|$. In this example, the PMC is

$2\Phi_{\beta}-\Phi_{\xi}=0.$ (8)

Under this condition, the QFI can be expressed as

$F_{\mathrm{m}}=2\bar{n}_{\mathrm{A}}\bar{n}_{\mathrm{B}}+\bar{n}_{\mathrm{A}}+\bar{n}_{\mathrm{B}}+2\bar{n}_{\mathrm{A}}\sqrt{\bar{n}_{\mathrm{B}}^{2}+\bar{n}_{\mathrm{B}}},$
(9)

where $\bar{n}_{\mathrm{A}}=|\beta|^{2}$ and
$\bar{n}_{\mathrm{B}}=\sinh^{2}|\xi|$. This equation is equivalent to the
corresponding equation in Ref. Jarzyna . In this expression of QFI, it is only
related to the average photon numbers in both ports: $\bar{n}_{\mathrm{A}}$
and $\bar{n}_{\mathrm{B}}$. Figure 2 shows the variation of
$F_{\mathrm{m}}/N^{2}$ with the change of $\bar{n}_{\mathrm{A}}$,
$\bar{n}_{\mathrm{B}}$. From this plot one can find that the optimal value of
quantum Fisher information for a fixed $N$ is obtained near the
$n_{\mathrm{A}}=n_{\mathrm{B}}$ line, especially for a large $N$. This is
because when $\bar{n}_{\mathrm{B}}$ is large, for fixed $N$,
$F_{\mathrm{m}}\simeq
4\bar{n}_{\mathrm{A}}\bar{n}_{\mathrm{B}}+\bar{n}_{\mathrm{A}}+\bar{n}_{\mathrm{B}}\leq
N^{2}+N$ and the optimal value can be achieved when
$\bar{n}_{\mathrm{A}}=\bar{n}_{\mathrm{B}}$. Adding the PMC, the optimal
choice for a large $N$ is that $\beta=\exp(i\Phi_{\xi}/2)\sinh|\xi|$. Also,
from this figure one can see that with the increase of $N$, the range that
$F_{\mathrm{m}}>N^{2}$ is increasing, which indicates that high intensity
input state is good for the enhancement of the phase sensitivity.

_PMC with unbanlanced beam splitter._ –For a more general beam splitter
transformation $\exp(i\tau J_{x})$ with $\tau\in[0,2\pi)$, the total setup of
the interferometer can be described by

$U_{\mathrm{MZ}}=\exp[i\theta(J_{z}\cos\tau-J_{y}\sin\tau)].$ (10)

If we restrict the input state $\rho_{\mathrm{A}}$ as a pure state, i.e.,
$|\psi\rangle\langle\psi|$, and $\rho_{\mathrm{B}}$ is still an even(odd)
state, the QFI reads $F=4\mathrm{Var}(J_{z}\cos\tau-J_{y}\sin\tau)$. Due to
the property of even(odd) state, one can find that

$F=4\cos^{2}\tau\mathrm{Var}(J_{z})+4\sin^{2}\tau\mathrm{Var}(J_{y}).$ (11)

It is known that for product states

$\mathrm{Var}(J_{z})=\frac{1}{4}\left[\mathrm{Var}(a^{\dagger}a)+\mathrm{Var}(b^{\dagger}b)\right],$
(12)

and $\mathrm{Var}(a^{\dagger}a)$, $\mathrm{Var}(b^{\dagger}b)$ are irrelevant
to the phase, therefore, the phase matching condition is decided by the second
item $\mathrm{Var}(J_{y})$ and keeps the same form with Eq. (4).

_PMC with photon losses._ –Now we consider the effects of photon losses on the
PMC. The scenario we use is shown in Fig. 1. Traditionally, the photon losses
can be described by a fictitious beam splitter Dobrzanski09 ; Dorner09 ;
Gkortsilas12 ; Huver08 ; Lee09 ; Xiang-Bin Wang . Here we use
$B^{T}=\exp[i(2\arccos\sqrt{T})J_{x}]$ to describe this fictitious beam
splitter. $T$ is the so-called transmission coefficient and we also define
$R=1-T$ as the reflection coefficient. When $T=1(R=0)$, there are no photon
losses in the interferometer and when $T=0(R=1)$, all the photons leak out of
the interferometer. This total photon loss scenario can be mapped into a neat
scenario, which includes two steps: first the input state goes through a
particle loss channel, then the output state imports into a MZ interferometer
without losses.

We assume that the leak in both arms share the same transmission coefficient
$T$ and the input state is separable:
$\rho_{\mathrm{A}}\otimes\rho_{\mathrm{B}}$. Then after the particle loss
channel, the reduced density matrix $\rho$ reads

$\rho=\mathrm{Tr}_{\mathrm{CD}}\left(\Gamma^{\mathrm{D}}\Gamma^{\mathrm{C}}\rho_{\mathrm{A}}\otimes\rho_{\mathrm{B}}\Gamma^{\mathrm{C}\dagger}\Gamma^{\mathrm{D}\dagger}\right),$
(13)

where the operator of the particle loss channel can be expressed by
$\Gamma^{\mathrm{C}}=\exp\left[i\sqrt{2}\arccos\sqrt{T}\left(J_{x}^{\mathrm{AC}}+J_{y}^{\mathrm{BC}}\right)\right]$
and
$\Gamma^{\mathrm{D}}=\exp\left[i\sqrt{2}\arccos\sqrt{T}\left(J_{x}^{\mathrm{BD}}+J_{y}^{\mathrm{AD}}\right)\right]$.
Here $J_{x}^{\mathrm{AC}}=\frac{1}{2}(a^{\dagger}c+ac^{\dagger})$ and
$J_{y}^{\mathrm{BC}}=\frac{1}{2i}(b^{\dagger}c-bc^{\dagger})$ with $c$,
$c^{\dagger}$ the annihilation and creation operators of mode C. So as
$J_{x}^{\mathrm{BD}}$ and $J_{y}^{\mathrm{AD}}$. Next $\rho$ goes through a
usual MZ interferometer, which can be described by $\exp(-\theta J_{y})$.
Utilizing this representation, this photon losses scenario can be classified
into a usual MZ interferometer scenario with a mixed input state.

Now we choose the initial states as $\rho_{\mathrm{A}}=|i\alpha
e^{i\Phi}\rangle\langle i\alpha e^{i\Phi}|$ and
$\rho_{\mathrm{B}}=|{\alpha}\rangle_{+}\langle{\alpha}|$. Here $\Phi$ is the
relative phase and $\Phi\in[0,\pi)$. Based on the PMC, the optimal value of it
is zero for no loss scenario. The similar scenario has been considered with
only one phase shift in one arm Joo ; G R Jin ; Jing . Through some
calculations, one can obtain the analytic expression of the QFI

$\displaystyle F=$ $\displaystyle
4T|\alpha|^{2}\\!\left[N_{\alpha}^{2}\\!+\\!T|\alpha|^{2}\\!\left(2N_{\alpha}^{2}-1\right)\\!\right]\\!+\\!4T^{2}|\alpha|^{4}\mathcal{G}\cos^{2}\Phi$
$\displaystyle-16T^{2}N_{\alpha}^{4}|\alpha|^{4}(1-p_{r}^{2})p_{t}^{2}\sin^{2}\Phi,$
(14)

where $p_{r}=\exp(-2|\alpha|^{2}R)$, $p_{t}=\exp(-2|\alpha|^{2}T)$ and
$\mathcal{G}=1-4N_{\alpha}^{4}\left(1-p_{r}^{2}\right)$. From this equation,
one can see that for a fixed $\alpha$ and transmission coefficient $T$, the
maximum value of Eq. (14) can always be obtained at $\Phi=0$, which indicates
that the phase matching condition remains _unchanged_ in this photon losses
case and is not effected by the transmission coefficient $T$. Under the phase
matching condition, the quantum Fisher information reads

$F_{\mathrm{m}}=4TN_{\alpha}^{2}|\alpha|^{2}+8T^{2}N_{\alpha}^{2}|\alpha|^{4}\left[1-2N_{\alpha}^{2}\left(1-p_{r}^{2}\right)\right].$
(15)

Utilizing the input average photon number $\bar{n}_{\mathrm{A}}$ and the input
total photon number $N$, it can be rewritten as

$F_{\mathrm{m}}=TN+2T^{2}N\bar{n}_{\mathrm{A}}\left[1-2N_{\alpha}^{2}\left(1-p_{r}^{2}\right)\right].$
(16)

This equation can reduce to Eq. (7) for $T=1(R=0)$. From this equation, one
can find that the photon losses have a negative influence on the QFI. With the
decrease of the transmission coefficient $T$, the QFI reduces. For a small
losses, namely $R$ is very small, the QFI reduces to

$F_{\mathrm{m}}=N+2N\bar{n}_{\mathrm{A}}-\left[1+4\bar{n}_{\mathrm{A}}(N+1)\right]NR.$
(17)

When $R<R_{c}$, with
$R_{c}=2\bar{n}_{\mathrm{A}}/\left[1+4\bar{n}_{\mathrm{A}}(N+1)\right]$, the
QFI can still larger than $N$, which indicates that these input states can
still surpass the shot-noise limit. This device is loss-tolerant and robust
within the region that $R\in[0,R_{c}]$. For a large $N$, $R_{c}\simeq 1/(2N)$.

_Conclusion_.–In summary, we have considered a general scenario of a Mach-
Zehnder interferometer, where the input state in one port of the
interferometer is an arbitrary state and the other one is an even(odd) state.
We have provided an analytic expression of the QFI and showed a general PMC
under which the phase sensitivity can be enhanced. We give two explicit
examples with one port in a coherent state and another in a even coherent
state or the squeezed vacuum state. The PMC plays an important role to improve
the sensitivity in these cases.

For the unbalanced beam splitters, we find the PMC is unchanged with one port
in a arbitrary pure stare and another in an even(odd) state. We also
considered the example of coherent superposition state with photon losses and
find that the PMC keeps unchanged and is not affected by the transmission
coefficients. Besides, for a small loss, we show that this setup is robust and
loss-tolerant. The present work sheds new light on the problem of how to
enhance phase sensitivity in quantum metrology.

_Acknowledgments_.–This work was supported by the NFRPC through Grant No.
2012CB921602, the NSFC through Grants No. 11025527 and No. 10935010.

## References

  * (1) V. Giovannetti, S. Lloyd, and L. Maccone, Nat. Photonics 5, 222 (2011); V. Giovannetti, S. Lloyd and L. Maccone, Phys. Rev. Lett. 96, 010401 (2006); C. C. Gerry and J. Mimih, Contemporary Physics, 51, 497-511 (2010).
  * (2) M. G. Genoni, S. Olivares, and M. G. A. Paris, Phys. Rev. Lett. 106, 153603 (2011); L. Pezzé and A. Smerzi, Phys. Rev. Lett. 100, 073601 (2008); Y. C. Liu, G. R. Jin and L. You, Phys. Rev. A 82, 045601 (2010); F. Benatti, R. Floreanini and U. Marzolino, J. Phys. B: At. Mol. Opt. Phys. 44, 091001 (2011); C. Gagatsos, O. Oreshkov and N. Cerf, Phys. Rev. A 87, 042307 (2013).
  * (3) K. P. Seshadreesan, P. M. Anisimov, H. Lee, and J. P. Dowling, New J. Phys. 13, 083026 (2011); D. W. Berry, B. L. Higgins, S. D. Bartlett, M. W. Mitchell, G. J. Pryde, and H. M. Wiseman, Phys. Rev. A 80, 052114 (2009); R. Birrittella, J. Mimih and C. C. Gerry, Phys. Rev. A 86, 063828 (2012).
  * (4) C. C. Gerry and A. Benmoussa, Phys. Rev. A 65, 033822 (2002); C. C. Gerry and R. A. Campos, Phys. Rev. A 64, 063814 (2001); B. M. Escher, R. L. de MatosFilho, and L. Davidovich, Nat. Phys. 7, 406 (2011);
  * (5) X.-M. Lu, S. Luo and C. H. Oh, Phys. Rev. A 86, 022342 (2012); X.-M. Lu, Z. Sun, X. Wang, S. Luo, and C. H. Oh, Phys. Rev. A 87 050302 (2013); G. Tóth and D. Petz, Phys. Rev. A 87, 032324 (2013).
  * (6) J. Estève, C. Gross, A. Weller, S. Giovanazzi, and M. K. Oberthaler, Nature (London) 455, 1216 (2008); M. F. Riedel, P. Böhi, Y. Li, T. W. Hänsch, A. Sinatra, and P. Treutlein, Nature (London) 464, 1170 (2010).
  * (7) J. Appel, P. J. Windpassinger, D. Oblak, U. B. Hoff, N. Kjaegaard, and E. S. Polzik, Proc. Natl. Acad. Sci. U.S.A. 106, 10960 (2009).
  * (8) T. Nagata, R. Okamoto, J. L. O’Brien, K. Sasaki, and S. Takeuchi, Science 316, 726 (2007).
  * (9) A. Fujiwara, Phys. Rev. A 63, 042304 (2001); A. Fujiwara, Phys. Rev. A 65, 012316 (2001); H. Imaiand A. Fujiwara, J. Phys. A 40, 4391 (2007); A. Fujiwara and H. Imai, J. Phys. A 36, 8093 (2003).
  * (10) A. Monras and M. G. A. Paris, Phys. Rev. Lett. 98, 160401 (2007); D. Collins, Phys. Rev A 87, 032301 (2013); J. Kahn, Phys. Rev. A 75, 022326 (2007); M. Hotta and T. Karasawa, Phys. Rev A 78, 012332 (2008) and so on.
  * (11) M. Ahmadi, D. E. Bruschi, N. Friis, C. Sabn, G. Adesso and I. Fuentes, e-print: arXiv:1307.7082 [quant-ph].
  * (12) C. M. Caves, Phys. Rev. D 23, 1693–1708 (1981).
  * (13) J. Joo, W. J. Munro and T. P. Spiller, Phys. Rev. Lett. 107, 083601 (2011).
  * (14) P. M. Anisimov, G. M. Raterman, A. Chiruvelli, W. N. Plick, S. D. Huver, H. Lee and J. P. Dowling, Phys. Rev. Lett. 104, 103602 (2010).
  * (15) L. Pezzé and A. Smerzi, Phys. Rev. Lett. 110, 163604 (2013).
  * (16) C. W. Helstrom, _Quantum Detection and Estimation Theory_ (Academic, New York, 1976).
  * (17) A. S. Holevo, _Probabilistic and Statistical Aspects of Quantum Theory_ (NorthHolland, Amsterdam, 1982).
  * (18) J. Liu, X. Jing, W. Zhong and X. Wang, submitted.
  * (19) B. Yurke, S. L. McCall and J. R. Klauder, Phys. Rev. A 33, 4033 (1986).
  * (20) R. A. Campos, B. E. A. Saleh and M. C. Teich, Phys. Rev. A 40, 1371 (1989).
  * (21) B.C. Sanders and G. J. Milburn, Phys. Rev. Lett. 75, 2944–2947 (1995).
  * (22) M. O. Scully and M. S. Zubairy, _Quantum Optics_ (Cambridge University Press, Cambridge, 1997).
  * (23) M. Jarzyna and R. Demkowicz-Dobrzański, Phys. Rev. A 85, 011801(R) (2012).
  * (24) R. Demkowicz-Dobrzanski, U. Dorner, B. J. Smith, J. S. Lundeen, W. Wasilewski, K. Banaszek, and I. A. Walmsley, Phys. Rev. A 80, 013825 (2009).
  * (25) S. D. Huver, C. F. Wildfeuer, and J. P. Dowling, Phys. Rev. A 78, 063828 (2008).
  * (26) U. Dorner, R. Demkowicz-Dobrzanski, B. J. Smith, J. S. Lundeen, W. Wasilewski, K. Banaszek, and I. A. Walmsley, Phys. Rev. Lett. 102, 040403 (2009).
  * (27) T. W. Lee, S. D. Huver, H. Lee, L. Kaplan, S. B. McCracken, C. Min, D. B. Uskov, C. F. Wildfeuer, G. Veronis, and J. P. Dowling, Phys. Rev. A 80, 063803 (2009).
  * (28) N. Gkortsilas, J. J. Cooper, and J. A. Dunningham, Phys. Rev. A 85, 063827 (2012).
  * (29) X.-X. Zhang, Y.-X. Yang, and X.-B. Wang, Phys. Rev. A 88, 013838 (2013).
  * (30) Y. M. Zhang, X. W. Li, W. Yang and G. R. Jin, e-print: arXiv:1307.7353.
  * (31) X. Jing, J. Liu, W. Zhong and X. Wang, e-print: arXiv:1307.8009.


# LongBet: Heterogeneous Treatment Effect Estimation in Panel Data

Meijia Wang
Google
Ignacio Martinez
Google
P. Richard Hahn
Arizona State University

###### Abstract

This paper introduces a novel approach for estimating heterogeneous treatment
effects of binary treatment in panel data, particularly focusing on short
panel data with large cross-sectional data and observed confoundings. In
contrast to traditional literature in difference-in-differences method that
often relies on the parallel trend assumption, our proposed model does not
necessitate such an assumption. Instead, it leverages observed confoundings to
impute potential outcomes and identify treatment effects. The method presented
is a Bayesian semi-parametric approach based on the Bayesian causal forest
model, which is extended here to suit panel data settings. The approach offers
the advantage of the Bayesian approach to provides uncertainty quantification
on the estimates. Simulation studies demonstrate its performance with and
without the presence of parallel trend. Additionally, our proposed model
enables the estimation of conditional average treatment effects, a capability
that is rarely available in panel data settings.

Keywords: Tree, BART, Bayesian Causal Forest, Panel Data, Heterogeneity,
Causal Inference, Staggered adoption, XBART, XBCF

##  1 Introduction

Causal inference problems in panel data involve understanding and quantifying
the causal relationships between the response and the treatment over time. It
is widely studied in policy evaluation and intervention analysis. This allows
policymakers to study the effect of a policy on target variables at different
time points, considering the temporal ordering of events. Panel data can
provide valuable insights into dynamic processes of the causal effect and
heterogeneity on different units.

Here we focus on the study of heterogeneous time-varying treatment effects in
panel data. Time-varying treatment effects refer to the situation where the
effect of a treatment or intervention may change over time, potentially
exhibiting different patterns or magnitudes at different points in time. The
key to the problem is to effectively estimate the potential outcomes for the
treated groups and identify the treatment effect given observed treated
outcomes.

Several methods have been developed in the Economics literature to study the
average treatment effects over time with different settings. Fixed effects
methods (Callaway and Karami,, 2023; Imbens et al.,, 2021) model the response
variable by fixed effects for individuals and factor loadings over time and
use it to impute the unobserved potential outcomes for the treated. Relying on
the parallel trend assumption, difference-in-differences design (Callaway and
Sant’Anna,, 2021; Borusyak et al.,, 2023; Abadie,, 2005) imputes the potential
outcomes for the treated parallel to observed control units based on the pre-
treatment differences. (Abadie,, 2021) provides an overview of the synthetic
control method, which predicts the potential outcome of each treated unit with
an ensemble of observed control units. Recent works (Acemoglu et al.,, 2013;
Kreif et al.,, 2016; Abadie and L’Hour,, 2021; Dube and Zipperer,, 2015) have
applied the synthetic control approach to a large number of units. Interrupted
time series analysis (Bernal et al.,, 2016, 2018; Schaffer et al.,, 2021)
focuses on studying the effect of population-level interventions by analyzing
the time series data before and after the intervention. It involves modeling
the trend and autocorrelation in the outcome variable while considering the
timing of the intervention. Other methods include the matrix completion
approach (Athey et al.,, 2021).

While most literature focuses on studying average treatment effects over time
for a treatment adopted at a particular time, recent works have been developed
to accommodate staggered adaption and estimate heterogeneous treatment
effects. (Borusyak et al.,, 2023) accommodate the heterogeneity in treatment
effect by leveraging the fixed effect model. (Callaway and Sant’Anna,, 2021)
relaxes the parallel trend assumption to hold only after conditioning on
observed covariates. (Wooldridge,, 2021) proposed a two-way Mundlak (TWM)
regression that allows for heterogeneity in treatment trends. (Xu,, 2017)
proposes a generalized synthetic control method that relaxes the parallel
trend assumption and allows for multiple treatment units and varying treatment
by combining synthetic control with a linear fixed effects model but requires
observed covariates to be time-varying. Xu also proposed other methods (to be
cited) that require time-varying covariates. (Athey et al.,, 2021; Tanaka,,
2020) propose using the matrix completion method to impute potential outcomes
as missing values and in turn estimate the treatment effect.

Traditional methods in causal inference often rely on the parallel trend
assumption, assuming that the untreated and treated groups follow parallel
trends over time. However, in various real-world scenarios, this assumption
may not hold, leading to limitations in the applicability of traditional
methods. For instance, in e-commerce, different verticals may exhibit
seasonality, and user-level data can vary significantly, making it challenging
to satisfy the parallel trend assumption. Similarly, in patient analysis,
individual responses to treatments may diverge, rendering the parallel trend
assumption inappropriate. Nevertheless, these datasets often provide valuable
observed covariates that can effectively characterize individual behaviors and
affect treatment responses. As a result, we propose a new method that can
leverage these observed covariates to impute potential outcomes and estimate
treatment effects, offering a powerful and flexible approach that can address
these complex scenarios without rigorously relying on the parallel trend
assumption.

(Hill,, 2011) first advocates using Bayesian Additive Regression Trees (BART)
in causal inference. BART’s flexibility and robustness to hyperparameters make
it an appealing choice for handling complex treatment effect heterogeneity.
Building on this foundation, Bayesian Causal Forest (BCF) (Hahn et al.,, 2020)
models and regularizes the heterogeneous treatment effects separately from the
prognostic effect of control variables using two BART forests in cross-
sectional settings. Additionally, (Krantsevich et al.,, 2022) introduced the
XBCF model which inherits the modeling of BCF with significantly accelerated
algorithm. Inspired by the BART literature, we propose a novel approach for
treatment effect estimation in panel data, thereby offering a unique and
distinctive approach that sets our work apart from existing literature in the
field.

Here we consider a scenario where we observe panel data of the response to a
binary treatment. We observe pre-treatment confoundings that affect the
treatment adoption and the treatment effect, which are observed. We consider a
”staggered rollout” design where the units can be treated at different times
but once a unit is treated it remains treated afterward. The goal is to
estimate the heterogeneous treatment effect conditioning on observed
confoundings. As a by-product, we can also provide an estimation of the
average treatment effect and/or average treatment effect on the treated by
averaging over the conditional treatment effect.

Our proposed method offers several distinct advantages over traditional
approaches in estimating time-varying treatment effects. Unlike methods that
rely on the parallel trend assumption, our model does not impose such
restrictive assumptions, allowing for greater applicability in scenarios with
diverse and complex data patterns. Moreover, our model is the first to harness
the power of Bayesian Additive Regression Trees (BART) for capturing time-
varying treatment effects, providing a novel and versatile framework. As a
semi-parametric model, it offers flexibility in capturing treatment effect
heterogeneity, accommodating variations in data across different time periods.
The model also readily handles staggered adoption of treatments, making it
suitable for a wide range of real-world applications. Leveraging a Bayesian
approach, our method provides credible intervals for conditional treatment
effects, enhancing the interpretability and reliability of the results.
Collectively, these features make our approach a powerful and innovative tool
for causal inference in the dynamic context of time-varying treatment effects.

##  2 Background

###  2.1 Bayesian additive regression trees (BART)

Bayesian Additive Regression Trees (BART), proposed by Chipman et al. (Chipman
et al.,, 2010), is a tree-based regression model known for its remarkable
capacity to model complex response surfaces, particularly in situations
characterized by large and intricate covariate structures. This model operates
through an ensemble of decision trees, providing both flexibility and
interpretability. In addition to its predictive power, BART is equipped with
the capability to offer Bayesian credible intervals, rendering it an
invaluable tool for conducting inference.

The BART model models the regression problem as the following:

$y_{i}=\sum_{l=1}^{L}g(x_{i};T_{l},\mu_{l})+\epsilon_{i},$ (1)

which seeks to approximate the relationship between covariates $\mathbf{x}$
and the response variable $y$ through a composition of regression trees. Here,
$T_{l}$ denotes a tree structure consisting of a set of partitioning rules
that divide the covariate space into distinct subspaces, represented as
$\\{\mathcal{A}_{1},\cdots,\mathcal{A}_{B}\\}$. Correspondingly,
$\mu_{l}=(\mu_{l1},\cdots,\mu_{lB})$ constitutes a vector of leaf parameters,
each associated with the leaf nodes within tree $T_{l}$.

$x_{1}\leq 0.8$$\mu_{1}$$x_{2}\leq 0.4$$\mu_{2}$$\mu_{3}$NoYesNoYes

0.40.8011$x_{1}$$x_{2}$$\mu_{1}$$\mu_{2}$$\mu_{3}$

Figure 1: A regression tree partitions a two-dimensional covariate space
(right) with two internal decision nodes (left). Each element of the partition
corresponds to a terminal node in the tree with leaf parameter $\mu_{l}$.

The binary regression tree shown in Figure 1 serves as an illustrative
example. This tree encompasses a collection of internal decision nodes as well
as terminal nodes, often referred to as leaf nodes. These leaf nodes
effectively partition the covariate space into a total of $B$ segments. The
function of the tree, denoted as $g(\mathbf{x},T_{l},\mu_{l})$, adopts a step-
like structure. This function allocates a point $\mathbf{x}$ to a leaf node
and assign the corresponding leaf parameter according to the split rules in
the tree $T_{l}$.

At the core of BART lies the formulation of its regularized tree prior and
leaf parameter prior, which characterizes its regularized nature while still
being flexible to model complex response surfaces. The tree prior put a strong
regularization on tree depth, the likelihood to split a node decreases
exponentially as the tree grows deeper. The leaf parameter prior assumes a
normal distribution, indicating that the model assumes that the responses
within each leaf node follow a Gaussian distribution. This encourages smoother
estimates within individual leaves, which, in turn, contributes to the overall
smoothness of the model’s predictions.

BART employs a random walk Metropolis-Hastings (MH-MCMC) algorithm to generate
trees from the posterior distribution. During each iteration, the algorithm
proposes a single growth or pruning procedure for each tree. The decision to
accept or reject these proposals is determined through a Metropolis-Hastings
(MH) approach, where the posterior probability of the proposed split is
evaluated. Furthermore, the iterative nature of BART involves fitting each
tree on the residuals of the ensemble formed by the other trees. This approach
ensures that subsequent trees focus on capturing the remaining patterns in the
data that have not been explained by the existing ensemble.

However, the algorithm’s efficiency diminishes when handling larger datasets
or covariate spaces. In contrast, a more recent approach known as XBART (He et
al.,, 2019; He and Hahn,, 2023) draws inspiration from BART but adopts a
different strategy. Instead of making incremental modifications to pre-
existing trees, XBART generates entirely new trees with each iteration. At
each node, the algorithm explores a comparable range of cutpoints as BART and
samples split criteria proportionate to the marginal likelihood. This
methodology facilitates the more efficient discovery of well-fitting deep
trees compared to BART MH-MCMC and performs better when handling larger
datasets. Our proposed method aligns with both XBART and BART frameworks,
although we will primarily demonstrate its efficacy using XBART for
computational efficiency.

###  2.2 Bayesian Causal Forest (BCF)

Constructed upon the foundation of BART model, the Bayesian Causal Forest
(BCF) model stands as a robust Bayesian non-parametric technique designed for
estimating heterogeneous treatment effects in cross-sectional data. Its
effectiveness has been demonstrated through applications on simulated and
real-world datasets (Hahn et al.,, 2020), showcasing superiority over popular
alternatives like BART, causal random forests, and propensity score matching
methods.

BCF identifies treatment effect based on the potential outcome framework
(Imbens and Rubin,, 2015) and relies on the assumption of the Stable Unit
Treatment Value Assumption (SUTVA) and the Strong Ignorability Assumption.
SUTVA stipulates that a unit’s potential outcome depends solely on its own
treatment status, unaffected by the treatments of others. Meanwhile, the
Strong Ignorability Assumption asserts that given observed covariates,
treatment assignment is independent of potential outcomes.

The BCF model formulates the response variable $Y$ as the sum of a prognostic
function and a treatment function given binary response $z_{i}$, represented
as:

$\mathbb{E}(Y_{i}|\mathrm{x}_{i},Z_{i}=z_{i})=\mu(\mathrm{x}_{i},\hat{\pi}(\mathrm{x}_{i}))+\tau(\mathrm{x}_{i})z_{i}.$
(2)

Here, $\mu$ captures the prognostic impact, and $\tau$ represents the
treatment effect. Both of these functions are estimated through BART forests,
which empowers the model to effectively approximate complex response surfaces.
Furthermore, (Hahn et al.,, 2020) suggested incorporating the estimated
propensity score $\hat{\tau}(\mathrm{x}_{i})$ in the prognostic effect
estimation, a strategy that has showcased enhanced accuracy.

XBCF (Krantsevich et al.,, 2022) is a recent innovation rooted in the Bayesian
Causal Forest (Hahn et al.,, 2020). This model expresses observable outcomes
as:

$\displaystyle y_{i}$
$\displaystyle=a\mu(\mathrm{x}_{i},\hat{\pi}_{i}))+b_{z_{i}}\tau(\mathrm{x}_{i})+\epsilon_{i},\quad\epsilon_{i}\sim
N(0,\sigma^{2}_{z_{i}}),$ (3) $\displaystyle a\sim N(0,1),\quad
b_{0},b_{1}\sim N(0,1/2)$

It extends the foundational Bayesian Causal Forest (BCF) technique with the
aim of boosting its performance in heterogeneous treatment effect estimation
and computational efficiency. While BCF employs the MH-MCMC algorithm from
BART to estimate both prognostic and treatment functions, XBCF leverages the
structure of XBART to enhance the efficiency of fitting the BCF model. The
model’s effectiveness is illustrated through comparisons with related
approaches in both simulation studies and empirical analyses.

The two models motivated us to extend their applicability to causal panel
data, thereby offering a unique and distinctive approach that sets our work
apart from existing literature in the field.

##  3 Methodology

In this section, we will begin by formalizing the problem, focusing on time-
varying treatment effect of a binary treatment in panel data. The binary
treatment can be applied to the treatment group at a specific time or be
implemented in a staggered roll-out design. Then, in order to identify
treatment effect, we will employt the dynamic potential outcomes framekwork
and extend the strong ignorability assumption to panel data. Finally, we will
introduce our proposed model LongBet.

Let $Y_{it}$ represent the continuous response of a unit $i$ at time $t$,
$Z_{it}$ indicate whether the unit $i$ has received the treatment at time $t$.
$\mathrm{X}_{i}$ represent a length $p$ vector of pre-treatment control
variables for the unit $i$. Let $T=\\{t_{0},\dots,t_{M}\\}$ denote the
observed period of time, $\mathrm{Y}_{i}$ and $\mathrm{Z}_{i}$ represent the
vector of responses and treatment indicators for the unit $i$, respectively.
We use lowercase letters, such as $y_{it}$ and vector $\mathrm{y}_{i}$, to
denote the observed values. The data will consist of $n$ independent
observations $(\mathrm{y}_{i},\mathrm{z}_{i},\mathrm{x}_{i})$.

In a staggered roll-out design, we consider a unit remains treated once it has
received treatment. Within this context, we define $S_{it}=\sum_{s\leq
t}Z_{is}$ to represent the time elapsed since treatment adoption, also
referred to as the post-treatment period in this study. Treatments assigned at
the same time would be a special case in the staggered roll-out design.

The dynamic potential outcomes framework was first introduced by (Robins,,
1986, 1987). This framework has been adapted for staggered adoption designs in
works such as (Sianesi,, 2004; Heckman et al.,, 2016; Callaway and Sant’Anna,,
2021). Under this framework, we use $Y_{it}(S)$ to denote the potential
outcome of unit $i$ at time $t$ after being treated for $S$ periods. In the
case where $S=0$, $Y_{it}(0)$ denotes the potential outcome without treatment.
The realized treatment effect is then calculated as

$Y_{it}=Y_{it}(0)+\sum_{S=1}^{T}\mathbbm{1}(S_{it}=S)(Y_{it}(S)-Y_{it}(0)).$
(4)

For our proposed method, we rely on the SUTVA assumption and extend the strong
ignorability assumption to the panel data settings. Namely, when there are no
unmeasured confounders, the potential outcomes are independent of the
treatment variable conditioning on observed confounders:

$Y_{it}(S)\perp\\!\\!\\!\perp\mathrm{Z}_{i}|\mathrm{X}_{i},\quad S=0,1,\dots,$
(5)

and that every individual has a non-zero probability of adopting the treatment
if it is not treated yet:

$0<\text{Pr}(Z_{it}=1|\mathrm{X}_{i},Z_{i(t-1)}=0)<1$ (6)

Following the above assumptions, the treatment effect originally defined by
the potential outcome framework can be relaxed to the difference in observable
outcomes conditioning on the confoundings and time since treatment:

$\displaystyle\tau_{t}(\mathrm{X}_{i},S)$
$\displaystyle:=\mathbb{E}[Y_{it}(S)|\mathrm{X}_{i}]-\mathbb{E}[Y_{it}(0)|\mathrm{X}_{i}]$
$\displaystyle=\mathbb{E}[Y_{it}|\mathrm{X}_{i},S_{it}=S]-\mathbb{E}[Y_{it}|\mathrm{X}_{i},S_{it}=0].$

Expanding from the XBCF model, we assume that the time-varying response to a
binary treatment can be approximated by the following equation:

$Y_{it}=\alpha\mu(\mathrm{X}_{i},T=t)+\beta_{S}\nu(\mathrm{X}_{i},S_{it},T=t)+\epsilon_{it}.$
(7)

The first term captures the prognostic effect in the absence of the treatment.
The second term models the dynamic treatment effect as a function of the
confoundings, the time since treatment and time $T$. We also introduce a
Gaussian process factor $\beta_{S}$ that captures the general trend of the
treatment effect shared among all units and varies by the time since
treatment. For simplicity, we assume the error term to be independent Gaussian
across units and time. The prognostic effect function $\mu(x_{i},t)$ and
treatment effect function $\nu(x_{i},t,s)$ are modeled by two separate XBART
forests, which are modified to handle panel data and consider splits on the
time dimensions.

Under the model specification, the treatment effect estimator is given by

$\tau_{t}(\mathrm{X}_{i},S)=\beta_{S}\nu(\mathrm{X}_{i},S_{it},T=t)-\beta_{0}\nu(\mathrm{X}_{i},S=0,T=t).$
(8)

The treatment effect estimate is obtained by averaging the posterior draws of
the estimator, and the method also provides credible intervals for the
treatment effect, conditioning on any covariates.

##  4 Simulation Study

To evaluate the performance of our proposed model, we adapt the simulation
studies from (Hahn et al.,, 2020) to generate panel data conditioning on
cross-sectional confoundings. We consider a total of $4$ data generating
process with the combination of $2$ types of prognostic effects (parallel and
non-parallel) and $2$ types of treatment effects (homogeneous and
heterogeneous).

We consider $5$ covariate variables which could affect the response variable
and treatment adoption but do not change over time. The first three covariate
variables are generated from the independent standard normal distribution
$x_{1},x_{2},x_{3}\overset{\text{i.i.d}}{\sim}N(0,1)$; the fourth is a binary
variable from binomial distribution; the last one is an unordered categorical
variable with three levels and equal probability.

The response variable consists of a prognostic effect, a treatment effect (if
being treated), and a noise term in the following form:

$y(\mathrm{x}_{i},t,s)=\eta(\mathrm{x}_{i},t)+\mathbbm{1}\\{s>0\\}\nu(\mathrm{x}_{i},t,s)+\epsilon_{it},\quad\epsilon_{it}\sim
N(0,0.25)$ (9)

When the parallel trend holds, the prognostic effect is the addition of a time
factor $f_{t}$ and a loading factor $\gamma(\mathrm{x})$ determined by cross-
sectional covariates. When the parallel trend is violated, the two factors
intertwine to determine the prognostic effect dynamically.

$\eta(\mathrm{x},t)=\begin{cases}f_{t}+\gamma(\mathrm{x})\quad&\text{parallel}\\\
f_{t}\gamma(\mathrm{x})\quad&\text{non-parallel.}\end{cases}$

The loading factor is defined as
$\gamma(\mathrm{x})=g(x_{4})++x_{1}|x_{3}-1|$, where $g(0)=2$ and $g(1)=-1$.
The time factor $f_{t}$ is generated from an ARIMA$(1,0,1)$ model with an
autoregressive coefficient (ar) of $0.7$ and a moving average coefficient (ma)
of $-0.4$ and $1$ is added to each observation to center the time factor
around the value 1.

The treatment effect is composed of a cohort effect, a time factor depending
on post-treatment time, and the covariates when considering heterogeneous
treatment effect.

$\nu(\mathrm{x},t,s)=\begin{cases}2c_{e_{i}}h_{s}\quad&\text{homogeneous}\\\
(2+x_{2}x_{5})c_{e_{i}}h_{s}\quad&\text{heterogeneous.}\end{cases}$

The cohort effect $c_{e_{i}}$ is determined by event time $e_{i}$ which is
when the unit $i$ adopts the treatment. The values of $c_{e_{i}}$ are listed
in Table 1 simulating the decaying treatment effect in many product adoption
scenarios. The time factor for the treatment effect is defined as
$h_{s}=se^{-s}$ to mimic the decaying treatment effect after adoption.

Event Time $e_{i}$ | 7 | 8 | 9 | 10 | 11 | 12
---|---|---|---|---|---|---
Cohort Effect $c_{e_{i}}$ | 1.3 | 1.2 | 1.1 | 1.0 | 0.9 | 0.8
Table 1: Cohort effect $c_{e_{i}}$ for each cohort that adopts the treatment
at time $e_{i}=7,8,\dots,12$.

Our simulations generate a total of $T=12$ time periods with $6$ pre-treatment
periods and the treatment starts rolling-out at $t=7$. Each untreated unit has
a chance of adopting the treatment based on their propensity score, which is
defined as

$\pi(\mathrm{x})=\frac{1}{5}\Phi(\frac{1}{2}\gamma(x)-\frac{1}{2}x_{1})^{2}+u_{i}/10,$

where $u_{i}\sim\text{Uniform}(0,1)$ and $\Phi(\cdot)$ is the quantile
function of standard normal. Starting at $t=7$, each untreated unit has the
probability of $\pi(\mathrm{x})$ of adopting the treatment.

We assess model performance from three aspects: 1) model accuracy; 2) model
reliability: how often the credible interval covers the truth; and 3) decision
impact: how well the model can inform decision making. Model accuracy is
evaluated through the average root mean squared error (RMSE) calculated across
ATT for each cohort and time point. We then assess the average coverage of the
95% simultaneous interval on the ATTs. In terms of decision impact, we examine
how frequently the confidence or credible intervals encompass zero for a non-
zero treatment effect, which serves as an indicator of the method’s capacity
to provide decision-relevant insights. For the LongBet model, we can further
measure the RMSE and coverage results specifically for the conditional average
treatment effect on the treated (CATT).111Under the strong ignorability
assumption, Longbet can estimate the average treatment effect. However, for
consistency with existing literature, we report only the estimated treatment
effects on the treated.

Let $ATT_{e,t}$ represent the average treatment effect on the treated (ATT)
for cohort $e$ at time $t$. This results in a total of $21$ ATT estimands
across $6$ cohorts and $6$ post-treatment periods. Denote the estimated ATT
for cohort $e$ at time $t$ as $\widehat{ATT}_{e,t}$. Furthermore, let
$q^{-}_{C,\alpha/2}\\{\widehat{ATT}_{e,t}\\}$ and
$q^{+}_{C,\alpha/2}\\{\widehat{ATT}_{e,t}\\}$ represents the lower and upper
bounds of the 95% simultaneous confidence interval (or credible interval for
Bayesian models) for the $C=21$ ATT estimands across cohorts and times. With
21 estimands in our scenario, the 95% simultaneous interval equates to a
99.76% interval for each estimand, as adjusted by Sidak. The evaluation
metrics for ATTs are formulated as follows. We report the metrics averaging
over 100 Monte Carlo simulations, offering a thorough analysis of the method’s
efficacy.

RMSE
$\displaystyle=\sqrt{\frac{1}{C}\sum_{e=t_{0}+1}^{t_{1}}\sum_{t=e}^{t_{1}}(ATT_{e,t}-\widehat{ATT}_{e,t})^{2}}$
Coverage
$\displaystyle=\frac{1}{C}\sum_{e=t_{0}+1}^{t_{1}}\sum_{t=e}^{t_{1}}\mathbbm{1}\left(q^{-}_{C,\alpha/2}\\{\widehat{ATT}_{e,t}\\}\leq
ATT_{e,t}\leq q^{+}_{C,\alpha/2}\\{\widehat{ATT}_{e,t}\\}\right)$ Cover0
$\displaystyle=\frac{1}{C}\sum_{e=t_{0}+1}^{t_{1}}\sum_{t=e}^{t_{1}}\mathbbm{1}\left(q^{-}_{C,\alpha/2}\\{\widehat{ATT}_{e,t}\\}\leq
0\leq q^{+}_{C,\alpha/2}\\{\widehat{ATT}_{e,t}\\}\right)$ $\displaystyle C$
$\displaystyle=\frac{2}{(t_{1}-t_{0}+1)(t_{1}-t_{0})}$

For LongBet, we set the number of iterations to $120$ (num_sweeps) and ${\tt
num\\_burnin}=20$. In each iteration, we train $50$ prognostic trees
(num_trees_pr) and $50$ treatment trees (num_trees_trt). We compare LongBet to
the following baseline methods:

  * •

DiD This is the difference-in-differences method proposed by (Callaway and
Sant’Anna,, 2021). It is tailored for staggered adoption design and relaxes
the parallel trend assumption to hold only after conditioning on observed
covariates. We employ their method’s doubly-robust estimand in our
simulations.

  * •

Non-linear DiD Considering the simulated prognostic function exhibits non-
linearity, we incorporate the known relationship $x_{1}x_{3}$ into the DiD
estimator to compare its performance with the original version.

  * •

DiD Imputation Introduced by (Borusyak et al.,, 2023), the DiD Imputation
method involves imputing potential outcomes using observed untreated responses
within a fixed effect model framework. This technique is capable of estimating
a weighted Average Treatment Effect on the Treated in scenarios featuring
heterogeneous causal effects, assuming parallel trends and the absence of
anticipation effects. 222The did_imputation package in R does not support
fitting with continuous cross-sectional covariates. Therefore, we proceeded to
fit the model solely with the inclusion of unit- and time-fixed effects..

Method | ATT | CATT
---|---|---
RMSE | Coverage | Cover 0 | RMSE | Coverage | Cover 0
Parallel, Homogeneous
LongBet | 0.022 | 0.944 | 0.078 | 0.031 | 0.963 | 0.149
DiD | 0.024 | 1.000 | 0.130 | - | - | -
Non-linear DiD | 0.024 | 0.999 | 0.132 | - | - | -
DiD Imputation | 0.017 | 0.999 | 0.095 | - | - | -
Parallel, Heterogeneous
LongBet | 0.019 | 0.965 | 0.074 | 0.068 | 0.979 | 0.308
DiD | 0.023 | 1.000 | 0.140 | - | - | -
Non-linear DiD | 0.023 | 1.000 | 0.140 | - | - | -
DiD Imputation | 0.017 | 0.999 | 0.099 | - | - | -
Non-parallel, Homogeneous
LongBet | 0.027 | 0.953 | 0.114 | 0.055 | 0.965 | 0.202
DiD | 0.314 | 0.616 | 0.228 | - | - | -
Non-linear DiD | 0.105 | 0.977 | 0.255 | - | - | -
DiD Imputation | 1.337 | 0.041 | 0.138 | - | - | -
Non-parallel, Heterogeneous
LongBet | 0.024 | 0.975 | 0.105 | 0.093 | 0.982 | 0.361
DiD | 0.296 | 0.690 | 0.239 | - | - | -
Non-linear DiD | 0.098 | 0.990 | 0.265 | - | - | -
DiD Imputation | 1.322 | 0.100 | 0.149 | - | - | -
Table 2: Simulation results

The simulation results, as summarized in Table 2, showcase the performance of
the examined methods in simulated scenarios. LongBet stands out as the
frontrunner in terms of RMSE on ATT. Remarkably, it achieves competitive
coverage on ATT while maintaining a low probability of covering 0. This
implies that LongBet is more likely to provide valuable insights for decision-
making purposes. Additionally, LongBet boasts low RMSE and high coverage on
CATT, reinforcing its proficiency. It’s worth noting that LongBet occasionally
presents lower coverage while achieving smaller RMSE, owing to its tighter
credible intervals that sometimes narrowly miss the true ATT.

In scenarios adhering to the parallel trend assumption, baseline methods
demonstrate competitive ATT estimation and higher coverage of the true
treatment effect. However, they also have higher coverage on $0$, indicating
lower statistical power to reject the “null hypothesis” - there is not enought
evidence to support the treatment effect is not zero.

On the other hand, when the assumption of parallel trends is not met, methods
such as DiD and DiD Imputation, which rely on this assumption, fail to
accurately estimate treatment effect. In contrast, the adapted version of DiD
demonstrates performance comparable to LongBet. However, there are certain
considerations to keep in mind. The non-linear DiD method requires prior
knowledge of the true relationship between the response variable and the
covariates, which is rarely available in real-world scenarios.

The simulation findings consistently emphasize LongBet’s effectiveness in
accurately estimating treatment effects and establishing credible intervals
across diverse scenarios. Moreover, these simulations bring to light LongBet’s
strength in situations where the prevalent parallel trend assumption, which
forms the basis of numerous existing methods for causal panel data, is not
applicable.

###  4.1 Prospective Assessment of Treatment Effects with LongBet

In the realm of causal inference, a comprehensive understanding extends beyond
historical treatment effects to include prospective analyses. A model’s
capability to anticipate future treatment effects can significantly influence
strategic planning and policy development. LongBet, our innovative method, is
designed to bridge this analytical divide by facilitating the forward
projection of treatment effects.

The LongBet framework encompasses various methodologies to project treatment
effects into forthcoming periods. The primary technique involves extrapolating
the treatment effect parameter $\beta$, employing Gaussian process regression
to discern temporal trends and patterns. Alternatively, LongBet can be
amalgamated with the local Gaussian process method detailed in (Wang et al.,,
2023), which utilizes treatment trees to enhance prediction accuracy. The
current discourse focuses on the primary extrapolation technique, evaluating
its predictive validity. Investigations into the latter method are earmarked
for subsequent research endeavors.

Figure 2: LongBet’s Average Treatment Effect Estimation and Forecasting. This
illustration delineates LongBet’s performance up to the 12th period (observed
data) and its subsequent forecasts for four additional periods. Each subplot
corresponds to a cohort characterized by its initial treatment period. The
black trajectory signifies the ground truth of treatment effects, while the
grey vertical dashed line demarcates the shift from historical data to
LongBet’s forecasts. The orange trajectory and the surrounding ribbon depict
LongBet’s treatment effect estimates and the accompanying $95\%$ credible
interval for each respective cohort.

This study augments the simulation setup previously introduced, extending the
observation window by four periods for future effect estimation. Each cohort’s
response to treatment is anticipated to exhibit a staggered effect, a
phenomenon deftly estimated by LongBet.

As visualized in Figure 2, the delineation between historical and forecasted
data is accentuated by a grey vertical line, providing a clear transition from
empirical data to LongBet’s extrapolations. The black line corresponds to the
ground truth. It is observed that during the observed periods, LongBet’s
predictions closely align with the ground truth, demonstrating its
effectiveness in delineating cohort-specific treatment trajectories.

In the predictive phase, LongBet’s estimates display an expansion of the
credible intervals, reflecting increased predictive uncertainty. However, for
cohorts treated later in the sequence, the intervals are notably narrower,
benefiting from the inference drawn from previous cohorts’ observations. We
have set the default parameters $\sigma=1$ and $\lambda=(t_{1}-t_{0})/2$ for
the Gaussian process extrapolation. These parameters can be fine-tuned,
potentially with domain expertise, to calibrate the prediction uncertainty.
Despite this, LongBet maintains estimations within a pragmatic range.

##  5 Conclusion

In this paper, we have introduced a novel approach for estimating the average
and conditional average treatment effects in the context of panel data. Unlike
conventional approaches in economics that heavily rely on the parallel trend
assumption, our proposed model does not necessitate such an assumption.
Instead, it leverages observed covariates to impute potential outcomes and
identify treatment effects. The method we presented is grounded in a Bayesian
semi-parametric framework, building upon the Bayesian causal forest model and
extended to suit the dynamics of panel data settings. The method combines the
strengths of Bayesian methodologies and offers credible intervals for
conditional average treatment effects.

It is important to acknowledge the underlying assumption of no unobserved
confounding in LongBet, which may not be universally applicable. However, the
proliferation of large-scale datasets in the era of big data offers potential
solutions to address this limitation. The semi-parametric nature of our
approach empowers it to capture complex and non-linear relationships between
covariates and responses, which can lead to more accurate treatment effect
estimation than traditional methods.

Looking ahead, our approach offers promising avenues for future research.
Extending LongBet to encompass time-varying covariates is a logical step
forward, allowing for a more nuanced understanding of how treatments interact
with changing covariate profiles over time. Additionally, extrapolating the
estimated causal effects by LongBet to provide predictive inferences for the
future could yield valuable insights for decision-makers.

In conclusion, LongBet emerges as a promising tool for causal inference in
panel data settings, offering a departure from conventional assumptions while
providing credible and insightful treatment effect estimates. As we
acknowledge its potential limitations and embrace its strengths, we envision a
dynamic landscape of research and application, where LongBet continually
evolves to address complex real-world scenarios and contributes to the
advancement of causal inference in panel data analysis.

## References

  * Abadie, (2005) Abadie, A. (2005).  Semiparametric Difference-in-Differences Estimators.  The Review of Economic Studies, 72(1):1–19.
  * Abadie, (2021) Abadie, A. (2021).  Using synthetic controls: Feasibility, data requirements, and methodological aspects.  Journal of Economic Literature.
  * Abadie and L’Hour, (2021) Abadie, A. and L’Hour, J. (2021).  A penalized synthetic control estimator for disaggregated data.  Journal of the American Statistical Association, 116(536):1817–1834.
  * Acemoglu et al., (2013) Acemoglu, D., Johnson, S., and Kermani, A. P. P. (2013).  The value of connections in turbulent times: Evidence from the united states.  Financial Crises eJournal.
  * Athey et al., (2021) Athey, S., Bayati, M., Doudchenko, N., Imbens, G., and Khosravi, K. (2021).  Matrix completion methods for causal panel data models.  Journal of the American Statistical Association, 116(536):1716–1730.
  * Bernal et al., (2016) Bernal, J. L., Cummins, S., and Gasparrini, A. (2016).  Interrupted time series regression for the evaluation of public health interventions: a tutorial.  International Journal of Epidemiology, 46:348 – 355.
  * Bernal et al., (2018) Bernal, J. L., Soumerai, S. B., and Gasparrini, A. (2018).  A methodological framework for model selection in interrupted time series studies.  Journal of clinical epidemiology, 103:82–91.
  * Borusyak et al., (2023) Borusyak, K., Jaravel, X., and Spiess, J. (2023).  Revisiting event study designs: Robust and efficient estimation.
  * Callaway and Karami, (2023) Callaway, B. and Karami, S. (2023).  Treatment effects in interactive fixed effects models with a small number of time periods.  Journal of Econometrics, 233(1):184–208.
  * Callaway and Sant’Anna, (2021) Callaway, B. and Sant’Anna, P. H. (2021).  Difference-in-differences with multiple time periods.  Journal of Econometrics, 225(2):200–230.  Themed Issue: Treatment Effect 1.
  * Chipman et al., (2010) Chipman, H. A., George, E. I., and McCulloch, R. E. (2010).  BART: Bayesian additive regression trees.  The Annals of Applied Statistics, 4(1):266 – 298.
  * Dube and Zipperer, (2015) Dube, A. and Zipperer, B. (2015).  Pooling multiple case studies using synthetic controls: An application to minimum wage policies.  IZA Discussion Paper, (8944).
  * Hahn et al., (2020) Hahn, P. R., Murray, J. S., and Carvalho, C. M. (2020).  Bayesian Regression Tree Models for Causal Inference: Regularization, Confounding, and Heterogeneous Effects (with Discussion).  Bayesian Analysis, 15(3):965 – 2020.
  * He and Hahn, (2023) He, J. and Hahn, P. R. (2023).  Stochastic tree ensembles for regularized nonlinear regression.  Journal of the American Statistical Association, 118(541):551–570.
  * He et al., (2019) He, J., Yalov, S., and Hahn, P. R. (2019).  Xbart: Accelerated bayesian additive regression trees.  In Chaudhuri, K. and Sugiyama, M., editors, Proceedings of the Twenty-Second International Conference on Artificial Intelligence and Statistics, volume 89 of Proceedings of Machine Learning Research, pages 1130–1138. PMLR.
  * Heckman et al., (2016) Heckman, J., Humphries, J., and Veramendi, G. (2016).  Dynamic treatment effects.  Journal of Econometrics, 191(2):276–292.
  * Hill, (2011) Hill, J. L. (2011).  Bayesian nonparametric modeling for causal inference.  Journal of Computational and Graphical Statistics, 20(1):217–240.
  * Imbens et al., (2021) Imbens, G., Kallus, N., and Mao, X. (2021).  Controlling for unmeasured confounding in panel data using minimal bridge functions: From two-way fixed effects to factor models.
  * Imbens and Rubin, (2015) Imbens, G. and Rubin, D. B. (2015).  Causal inference for statistics, social, and biomedical sciences: Sensitivity analysis and bounds.
  * Krantsevich et al., (2022) Krantsevich, N. M., He, J., and Hahn, P. R. (2022).  Stochastic tree ensembles for estimating heterogeneous effects.  In International Conference on Artificial Intelligence and Statistics.
  * Kreif et al., (2016) Kreif, N., Grieve, R., Hangartner, D., Turner, A. J., Nikolova, S., and Sutton, M. (2016).  Examination of the synthetic control method for evaluating health policies with multiple treated units.  Health Economics, 25(12):1514–1528.
  * Robins, (1986) Robins, J. M. (1986).  A new approach to causal inference in mortality studies with a sustained exposure period—application to control of the healthy worker survivor effect.  Mathematical Modelling, 7:1393–1512.
  * Robins, (1987) Robins, J. M. (1987).  Addendum to “a new approach to causal inference in mortality studies with a sustained exposure period—application to control of the healthy worker survivor effect”.  Computers & Mathematics With Applications, 14:923–945.
  * Schaffer et al., (2021) Schaffer, A. L., Dobbins, T. A., and Pearson, S.-A. (2021).  Interrupted time series analysis using autoregressive integrated moving average (arima) models: a guide for evaluating large-scale health interventions.  BMC Medical Research Methodology, 21.
  * Sianesi, (2004) Sianesi, B. (2004).  An evaluation of the swedish system of active labor market programs in the 1990s.  The Review of Economics and Statistics, 86(1):133–155.
  * Tanaka, (2020) Tanaka, M. (2020).  Bayesian matrix completion approach to causal inference with panel data.
  * Wang et al., (2023) Wang, M., He, J., and Hahn, P. R. (2023).  Local gaussian process extrapolation for bart models with applications to causal inference.  Journal of Computational and Graphical Statistics, 0(0):1–12.
  * Wooldridge, (2021) Wooldridge, J. (2021).  Two-way fixed effects, the two-way mundlak regression, and difference-in-differences estimators.  SSRN Electronic Journal.
  * Xu, (2017) Xu, Y. (2017).  Generalized synthetic control method: Causal inference with interactive fixed effects models.  Political Analysis, 25(1):57–76.