\section{Introduction}
Black holes are arguably the most fascinating predictions of Einstein's general theory of relativity. It is now evident that both astronomical observations and the theoretical formulation of general relativity suggest their real existence in nature. The most striking evidence for their existence is the detection of the
gravitational waves \cite{gw} and the shadow image of the supermassive M87*  \cite{eh1,eh2}. Black holes are extremely important in explaining many other problems in astrophysics that have to do with the high-energy phenomena in the form of X-ray emission and jets, accretion of matter, quasi-periodic oscillations, and the motions of nearby objects in orbit around the hidden mass \cite{GRAVITY:2020gka}.

Despite all these successes, there are still some conceptual problems related to black holes in general relativity. More specifically, the main problem associated with black holes is the existence of singularities at the coordinate center. As of today, this problem has not been solved. Many physicists have speculated that a quantum theory of gravity can solve this problem. However, the phenomenology of black holes remains a hot topic after the EHT findings. The black hole mimickers such as regular black holes, wormholes, naked singularities, grava-stars, etc. are also growing significantly in the literature \cite{Hayward:2005gi,Visser:2003ge,joshi,Gao:2016bin,Mosani:2020ena,Mosani:2020mro,Mosani:2021,az,amir}.

The concept of a regular black hole was first introduced by Bardeen in 1968 \cite{Bardeen:1968}, see also the more recent references \cite{Bergmann-Roman,Hayward:2005,Bardeen:2014,Frolov:2014,Frolov:2014,Frolov:2014b,Frolov:2016,Frolov:2017,Frolov:2018,Cano:2018,Bardeen:2018,regular,beyond}, and it has been intuitively interesting due to its non-singular structure.  One motivation for exploring such compact objects is that they could be black hole alternatives. Furthermore, as is widely known, many of the possible observational signatures of the black holes such as precession of timelike bound orbits, gravitational lensing and shadow properties can be mimicked by black hole mimickers \cite{Bambhaniya:2019pbr,Joshi:2019rdo,Bambhaniya:2020zno,Bambhaniya:2021ybs,Nampalliwar:2021tyz,Jusufi:2020zln,Saurabh:2020zqg,Nampalliwar:2021oqr,Joshi:2020tlq,Dey:2019fpv}.

From a historical perspective, it's worth noting that Einstein and Rosen proposed the existence of "bridges" through space-time \cite{ERb} using the theory of general relativity. These bridges connect two places in space-time, allowing for the creation of Einstein-Rosen bridges, also known as wormholes. Such wormholes have been proved to be non-traversable. The notion of having traversable wormholes was later investigated in the pioneering work of Morris and Throne \cite{Morris}.

The exotic matter that satisfies the flare-out criterion and violates the weak energy condition is required to sustain the structure of a wormhole (see \cite{Morris:1988cz, Morris:1988tu}). However, quantum matter fields have recently been proved to provide enough negative energy to allow some wormholes to be traversed. As a result, an exotic matter with negative energy density and large negative pressure, with a greater value than the energy density is required to form such a traversable wormhole \cite{k1,k2}.

We can therefore try to test black holes and wormholes using astrophysical data. In Ref. \cite{k3} authors studied the possibility to test wormholes geometries in our galaxy using infalling and radiating gas accretion model and the motion of the S2 star. The possibility of testing wormholes by observing the galactic center using the infalling gas model was studied in Refs. \cite{b1,b2}. Other interesting works concerning the formation and the stability of wormholes can be found in Refs. \cite{s1,s2,s3}. Based on astrophysical observations, today we know that in the galactic center of many galaxies a disk is formed. Namely, such a disk is made of rapidly rotating gas that slowly spirals onto a compact body which is assumed to be a black hole. During such accretion, the gravitational energy through the friction of the heat converts into radiation, which partially escapes, and cools down the
accretion disc. The only information that we have about
accretion disk physics comes from this radiation, when it
reaches radio, optical and X-ray telescopes, allowing astronomers to analyze its electromagnetic spectrum, and
its time variability. The radiation
properties of thin accretion disks were further analyzed
in the recent papers \cite{Shaikh:2019hbm,1,2,3,4,5,6}, where the effects of photon capture by the black hole on the spin evolution were presented as well.

Simpson and Visser proposed a very simple theoretically appealing spherically symmetric and static spacetime family \cite{Simpson:2018tsi}, derived from Schwarzschild geometry, that allows for a unique description of regular black hole and wormholes by smooth interpolation between these two possibilities using a length-scale parameter $'l'$ that drives the regularization of the central singularity. Its rotational form was recently shown as well \cite{Mazza:2021rgq}. In the present work, we are interested to use the Simpson and Visser metric and explore the NT accretion disk properties and images. It will be interesting to see whether or not we can distinguish the regular/non-regular black hole from the wormhole case.  This metric is of particular interest since it can describe a regular/non-regular black hole or traversable wormhole depending on the value of the regularization parameter $l$. Note that many physical properties of this metric, including the effect of charge, gravitational lensing, quasi-normal modes, and other effects were studied by many authors \cite{sv1,sv2,sv3,sv4,sv5,sv6,sv7,sv8}. In the present work, we are interested to study the optical appearance and the electromagnetic properties of the NT accretion disc around the Simpson-Visser metric.

This paper is organized as follows. In Section II, we present the Simpson-Visser metric. In Section III, we study the embedding diagram. In Section IV, we study the NT accretion disk and the images. In Section V, we explore the electromagnetic properties of the accretion disk. Finally, in Section VI, we comment on our results. Throughout this paper, we consider gravitational constant $(G)$ and speed of light $(c)$ as unity.

\section{Simpson-Visser spacetime}
The proposed spherically symmetric and static Simpson-Visser metric specifies a class of black hole mimickers with a minimal surface in place of the central singularity. The geometry of spacetime is given by \cite{Simpson:2018tsi},
\begin{equation}
ds^2=-f(r)dt^2+\frac{dr^2}{f(r)}+(r^2+l^2)(d\theta^2+\sin^2\theta d\phi^2).
\label{SVmetric}
\end{equation}
where,
\begin{equation}
    f(r)=\left(1-\frac{2M}{\sqrt{r^2+l^2}}\right).\; \; \; \;

\end{equation}
In the above spacetime, $M\ge 0$ represents the ADM mass, and $l>0$ is a parameter responsible for the regularisation of the central singularity and
possibly reflecting the quantum gravity effects. The above spacetime (\ref{SVmetric}) neatly interpolates between the standard Schwarzschild black hole and the Morris–Thorne traversable wormhole; Passing through a black-bounce (into a future incarnation of the universe), an extremal null-bounce (into a future incarnation of the universe), and a traversable wormhole are all intermediary stages. Therefore, this spacetime consist family of solutions as the black hole mimickers, which can be described as,
\begin{enumerate}
\itemsep-1pt
\item
The ordinary Schwarzschild spacetime ($l=0$);
\item
A ``black-bounce'' (regular black hole) with a one-way space-like throat  ($l<2M$);
\item
A one-way wormhole with a null throat ($l=2M$),
\item
A traversable wormhole in the Morris--Thorne sense ($l>2M$)
\end{enumerate}
The Simpson and Visser metric is everywhere regular as long as the parameter $l$ is non-zero, resulting in an unusual kind of "regular black hole," where the"origin" at r = 0 can be space-like, null, or timelike. Moreover, the Carter–Penrose diagrams and curvature tensors are defined in \cite{Simpson:2018tsi} and shown that the Einstein tensor has non-zero (mixed) components. In the present work we are interested to understand the effect of $l$ on the physical properties of the electromagnetic waves emitted by the accretion disk.

\begin{figure}
    \centering
    \includegraphics[scale=0.79]{emb1.png}
    \caption{The geometry of Simpson-Visser wormhole embedded in a three-dimensional Euclidean space. We have used $l=2.1$ and $M=1$.}
    \label{geometry}
\end{figure}

\begin{figure*}
    \centering
    \includegraphics[scale=0.69]{Veff1.pdf}
        \includegraphics[scale=0.69]{rms1.pdf}
    \caption{Left panel: The effective potential for different values of $l$ and $M=1$.
 Right panel: The radius of the marginally stable circular orbit or the ISCO for different values of $l$ and $M=1$. }
\end{figure*}

\section{Embedding Diagram}
	Let us study the geometry of the Simpson-Visser metric by embedding it into a higher-dimensional Euclidean space. To simplify the problem let us consider the  equatorial plane $\theta=\pi/2$ at a fixed moment  $t=$ Constant, in that case we have
	\begin{equation}
	ds^2=\frac{dr^2}{1-\frac{b(r)}{r}}+\mathcal{R}^2d\phi^2, \label{emb}
	\end{equation}
	where
	\begin{equation}
	b(r)= r\left(1-f(r)\right),\,\,\,\mathcal{R}^2=r^2+l^2.
	\end{equation}
	Let us embed this black hole metric into three-dimensional Euclidean space  in the cylindrical coordinates,
		\begin{eqnarray}
		ds^2&=&dz^2+d\mathcal{R}^2+\mathcal{R}^2d\phi^2,
		\end{eqnarray}
		which can be written as follows
		\begin{eqnarray}
		ds^2&=&\left[\left(\frac{d\mathcal{R}}{dr}\right)^2+\left(\frac{dz}{dr}\right)^2 \right]dr^2+\mathcal{R}^2d\phi^2.
		\end{eqnarray}
		From these equations,  we find that
	\begin{equation}
	\frac{dz}{dr}=\pm \sqrt{\frac{r}{r-b(r)}-(\frac{d\mathcal{R}}{dr})^2},
	\end{equation}

 As a particular case we will consider the traversable wormhole with $l>2$. Note that the integration of the last expression cannot be accomplished analytically. Invoking numerical techniques allows us to illustrate the embedding diagrams for the Simpson-Visser traversable wormhole given in figure (\ref{geometry}).

\begin{figure*}
    \centering
    \includegraphics[scale=0.50]{panel1.pdf}
    \caption{Images of the Novikov-Thorne thin-accretion disk around Schwarzschild Black hole and Simpson-Visser regular Black hole with $l=1$ and wormhole with $l=2.5$ for inclination 90$\degree$.}
    \label{fig:NT}
\end{figure*}

\begin{figure*}
    \centering
    \includegraphics[scale=0.50]{panel2.pdf}
    \caption{Images of the Novikov-Thorne thin-accretion disk around Schwarzschild Black hole and Simpson-Visser regular Black hole with $l=1$ and wormhole with $l=2.5$ for inclination 45$\degree$.}
    \label{fig:NT2}
\end{figure*}

\section{Novikon-Thorne thin accretion disks}
Massive particles move in stable circular timelike geodesics in the geometrically thin accretion disk. The Novikov-Thorne model of a thin accretion disk is considered here. There are two constants of motion along the timelike geodesics because we are dealing with spherically symmetric and static spacetimes. The particles' energy $(\Tilde{E})$ and angular momentum $(\Tilde{L})$ per unit rest mass, which are represented as,
\begin{equation}
    \Tilde{E}=\left(1-\frac{2M}{\sqrt{r^2+l^2}}\right)\Dot{t}, \; \; \; \; \; \; \;
    \Tilde{L}=(r^2+l^2)\Dot{\phi},
   \label{conserved}
\end{equation}
where over dot represents the derivative with respect to proper time $(\tau)$.
Using the above conserved quantities and timelike geodesic condition $(ds^2=-m^2)$, we can write the general form of spherically symmetric and static spacetime (\ref{SVmetric}) as,
\begin{equation}
    \Tilde{E}^2=\Dot{r}^2+ V_{eff}(r),
\end{equation}
this above equation represents the total energy of the particle, where $V_{eff}(r)$ is the effective potential of the Simpson-Visser spacetime for timelike geodesics. The effective potential of the Simpson-Visser spacetime is defined as,
\begin{equation}
    V_{eff}(r)=\Tilde{L}^2\frac{\left(1-\frac{2M}{\sqrt{r^2+l^2}}\right)}{(r^2+l^2)}+\left(1-\frac{2M}{\sqrt{r^2+l^2}}\right)m^2.
\end{equation}
where, for photons $m = 0$ and for test particles $m > 0$ (see Fig. 2). The sign of the circular geodesic orbits determines their stability. The condition $V''_{eff}(r_c)=0$ determines the marginally innermost stable circular orbit (ISCO). It indicates the ISCO position at the radius,
\begin{eqnarray}
r_{\rm ISCO}=\sqrt{(6M)^2-l^2}.
\label{isco}
\end{eqnarray}
The horizon position of a regular black hole $(l<2M)$ is defined as,
\begin{equation}
    r_h=\sqrt{(2M)^2-l^2}.
\end{equation}
For the shadow size, one can calculate the radius of the circular orbit of a photon $(l<3M, ds^2=0)$ as,
\begin{equation}
    r_{ph}=\sqrt{(3M)^2-l^2},
\end{equation}
and the critical impact parameter $b = E/L$ of the photon circular orbits is given by the expression,
\begin{equation}
    b^2_{ph}=27M^2,
\end{equation}
having the same size as a photon circular orbit around a Schwarzschild black hole and does not depend on $l$ \cite{sv7}. Note that for $l>3M$, there is no photon sphere but we have a contribution from the wormhole throat.
Now, the stable circular timelike geodesic can be obtained by satisfying the conditions,
\begin{equation}
    V_{eff}(r)=\Tilde{E}^2, \; \; \; \;
    V'_{eff}(r)=0, \; \; \; \;
     V''_{eff}(r)<0,
\end{equation}
where, the first two conditions are useful to get $\Tilde{E}$ and $\Tilde{L}$ as a function of radius of the circular orbit as,
\begin{equation}
    \Tilde{E}=\frac{1-\frac{2M}{\sqrt{r^2+l^2}}}{\sqrt{1-\frac{2M}{\sqrt{r^2+l^2}}-(r^2+l^2)\Omega^2}},
    \label{Etilta}
\end{equation}
\begin{equation}
    \Tilde{L}=\frac{(r^2+l^2)\sqrt{\frac{M}{(r^2+l^2)^{3/2}}}}{\sqrt{1-\frac{2M}{\sqrt{r^2+l^2}}-(r^2+l^2)\Omega^2}},
    \label{Ltilta}
\end{equation}
where, $\Omega$ gives the angular momentum of the particles in terms of the equation \cite{Shaikh:2019hbm},
\begin{equation}
    \Omega=\frac{d\phi}{dt}=\sqrt{\frac{f'(r)}{2r}}=\sqrt{\frac{M}{(r^2+l^2)^{3/2}}},
    \label{Omega}
\end{equation}
with prime representing a derivative with respect to the radial coordinate `r'.

\begin{figure}
    \centering
        \includegraphics[width=\columnwidth]{cross_sec_inten.png}
    \caption{The plot shows the horizontal cross-sectional intensity of the Novikov-Thorne accretion disks for $l=0$ and $l=1.25$, respectively.}
\end{figure}

\begin{figure*}
    \centering
        \includegraphics[scale=0.68]{Tem2.pdf}
    \includegraphics[scale=0.70]{Flux2.pdf}
    \caption{Left panel: The temperature profile of the thin accretion disk for different values of $l$ using the case $M=10^6 M_0$ and $\dot M_0 \sim 2.5\,\times 10^{-5} M_{\odot} /{\rm yr}$. Right panel: The  radiated energy flux over the thin accretion disk for different values of $l$ and $M=1$. In both plots we have $l=0$ (black curve), $l=1.8$ (red curve), $l=3.0$ (blue curve) and $l=4.0$ (orange curve).}
\end{figure*}

\begin{figure*}
    \centering
        \includegraphics[scale=0.65]{Tem1.pdf}
    \includegraphics[scale=0.72]{Flux1.pdf}
    \caption{Left panel: The temperature profile of the thin accretion disk for different values of $l$ using the case $M=10^6 M_0$ and $\dot M_0 \sim 8\,\times 10^{-12} M_{\odot} /{\rm yr}$. Right panel: The  radiated energy flux over the thin accretion disk for different values of $l$ and $M=1$. In both plots we have $l=0$ (black curve), $l=1.8$ (red curve), $l=3.0$ (blue curve) and $l=4.0$ (orange curve) in both plots.}
\end{figure*}

\begin{figure*}
    \centering
            \includegraphics[scale=0.69]{Lum22.pdf}
            \includegraphics[scale=0.71]{Lum.pdf}
    \caption{Left panel: The observed luminosity of the thin accretion disk for different values of $l$ using the case $M=10^6 M_0$ and $\dot M_0 \sim 2.5\,\times 10^{-5} M_{\odot} /{\rm yr}$.  Right panel: The observed luminosity of the thin accretion disk for different values of $l$ using the case $M=10^6 M_0$ and $\dot M_0 \sim 8\,\times 10^{-12} M_{\odot} /{\rm yr}$. We have $l=0$ (black curve), $l=1.8$ (red curve), $l=3.0$ (blue curve) and $l=4.0$ (orange curve) in both plots. }
\end{figure*}

In this part, we use the Novikov-Thorne model of a thin accretion disc composed of anisotropic fluid moving in the equatorial plane to examine the steady-state thin accretion disc. Certain structure equations govern the physical properties of the accretion disc, which arises from the necessity of the fluid's rest mass, energy, and angular momentum conservation. In this setup, the accreting matter in the disk can be described by the stress-energy tensor given by  \cite{1}.
\begin{eqnarray}
T^{\mu\nu}=\rho_0 u^\mu u^\nu +2 u^{(\mu}q^{\nu)}+t^{\mu\nu},
\end{eqnarray}
where
\begin{eqnarray}
u_\mu q^\mu&=&0,\\
u_\mu t^{\mu\nu}&=&0.
\end{eqnarray}
In the above equations we have the introduced the following quantities: $\rho_0$, $q^\mu$, and $t^{\mu\nu}$, which are known as the rest mass density of the accreting matter, the energy flow vector, and the stress tensor, respectively. These quantities are specified in the averaged rest-frame of the particle described by the 4-velocity $u^\mu$. If the rest mass is conserved, the equation can be used as,
\begin{equation}
\nabla_\mu (\rho_0 u^\mu)=0
\end{equation}
As a result, the disc radius has no effect on the time averaged rate of rest mass accretion.
\begin{eqnarray}
\dot{M}_0 &=& -2\pi \sqrt{-g} \Sigma u^r={\rm const.},
\end{eqnarray}
where $\Sigma(r)$ denotes the time-averaged surface density in cylindrical coordinates, and
\begin{eqnarray}
\Sigma (r) &=& \int^H_{-H}<\rho_0>dz,
\end{eqnarray}

The law of energy conservation and the law of angular momentum conservation both state that
\begin{eqnarray}
\nabla_\mu T^{t \mu}&=&0,\\
\nabla_\mu T^{\phi \mu}&=&0.
\end{eqnarray}

These relationships can be used to obtain the thin disk's time-averaged radial structure equations,
\begin{eqnarray}
&&[\dot{M}_0\tilde{E}-2\pi \sqrt{-g} \Omega W^r_\phi]_{,r} =4 \pi r F(r)\tilde{E}, \\
\end{eqnarray}
and
\begin{eqnarray}
&&~~[\dot{M}_0\tilde{L}-2\pi \sqrt{-g} W^r_\phi]_{,r}=4 \pi r F(r) \tilde{L}.
\end{eqnarray}
where $W^r_\phi$ is the averaged torque and is given by
\begin{eqnarray}
W^r_\phi=\int^H_{-H}<t^r_\phi>dz.
\end{eqnarray}
At this point, it's important to emphasise that the $\phi$-$r$ component of the stress tensor averaged over a characteristic time scale $\Delta t$  and the azimuthal angle $\Delta \phi=2 \pi$ is the quantity $<t^r_\phi>$. The flux $F(r)$ of the radiant energy over the disc can be represented in terms of the specific energy, angular momentum, and angular velocity of the orbiting particle in the thin accretion by using the energy-angular momentum relation $\tilde{E}_{,r}=\omega \tilde{L}_{,r}$. The flux of electromagnetic radiation emitted from a radial position 'r' of a disc is given by the standard formula, which may be deduced from the equations (\ref{Etilta}),(\ref{Ltilta}), and (\ref{Omega}) \cite{Shaikh:2019hbm},
\begin{equation}
    F(r)=-\frac{\Dot{M}}{4\pi\sqrt{-g}}\frac{\Omega'}{(\Tilde{E}-\Omega\Tilde{L})^2}\int_{r_{in}}^{r} (\Tilde{E}-\Omega\Tilde{L})\Tilde{L'} \,dr ,
\end{equation}
Note that $r_{in}$ represent the inner edge of the disk and $\Dot{M}$ is the mass accretion rate. In Figs. 3 and 4, we show the optical appearance of the Simpson-Visser spacetime surrounded by NT accretion disk. We use the ray-tracing formalism described in \cite{Younsi:2016azx} to numerically integrate the geodesic equation using \texttt{RK45} method with adaptive step size.
One can see that, the images are almost similar compared to the Schwarzschild geometry, which makes it difficult to distinguish these two spacetime geometries. However, in what follows we shall study in more details this problem, in particular we shall study in the electromagnetic properties of the Simpson-Visser spacetime and compere to the Schwarzschild geometry.

\section{The Radiating flux, temperature profile, and the luminosity of thin accretion disk}
In this section, we turn our attention to explore in more details the electromagnetic properties of the accretion disk, such as the radiation flux $F(r)$, the temperature profile, luminosity, and the accretion efficiency. Toward this purpose, we will numerically analyse the behaviour of the above quantities with respect to the parameter $l$. Furthermore we assume a central compact object with mass $M=10^6 M_{\odot}$ (a supermassive black hole), along two cases for the accretion rate given as follows \cite{1,6}
\begin{equation}
\dot M_0 \sim 10^{-5} M_{\odot} /{\rm yr}\,\, \text{and} \,\,\dot M_0 \sim 10^{-12} M_{\odot} /{\rm yr}
\end{equation}

Recall that the radiation flux $F(r)$ given by Eq. (4.24) represent the emitted flux by the thin accretion disk in the Simpson-Visser spacetime. Moreover, we expect that the Stefan-Boltzmann law holds and, therefore, in terms of the effective temperature associated to the accretion disk we have
\begin{eqnarray}
T_{\rm eff}(r) = \left(\frac{F(r)}{\sigma}\right)^{1/4},
\end{eqnarray}
where for the Stefan-Boltzmann constant we have $\sigma = 5.67 \times 10^{-5}\; {\rm erg}\; s^{-1} \;{\rm cm}^{-2}\; K^{-4}$. Our analyses shows that, the temperature and the radiated energy flux of the disk increases by increasing the parameter $l$, as can be seen from Figs. 6 and 7. As a specific example, take the case $\dot M_0 \sim 8\,\times 10^{-12} M_{\odot} /{\rm yr}$, we find that the maximal value for the temperature at a particular $l$ occurs at:
\begin{equation}\notag
l=0.0:\,T_{max}\simeq 867.01\, {\rm K}, \,\,r=8.70[M]
\end{equation}
\begin{equation}\notag
l=1.8:\,T_{max}\simeq 871.79 \, {\rm K},\,\,r=8.42 [M]
\end{equation}
\begin{equation}\notag
l=3.0:\,T_{max}\simeq 881.04\, {\rm K}\,\,r=8.09 [M]
\end{equation}

On the other hand, the maximal value for the radiating flux at particular $l$ occurs at:
\begin{equation}\notag
l=0.0:\,F_{max}\simeq 3.203 \times 10^7 {\rm erg}\; s^{-1} \;{\rm cm}^{-2}, \,r=8.70[M]
\end{equation}
\begin{equation}\notag
l=1.8:\,F_{max}\simeq 3.275 \times 10^7 {\rm erg}\; s^{-1} \;{\rm cm}^{-2},\,r=8.42 [M]
\end{equation}
\begin{equation}\notag
l=3.0:\,F_{max}\simeq 3.416 \times 10^7 {\rm erg}\; s^{-1} \;{\rm cm}^{-2},\,r=8.09 [M]
\end{equation}

The observed luminosity $L(\nu)$ of the thin accretion disc around the Simpson-Visser spacetime has a red-shifted black body spectrum while we assume the radiation emitted by the thin accretion disc surface to be perfect black body radiation \cite{1},
\begin{eqnarray}
L(\nu) &=& 4 \pi d^2 I(\nu)
\end{eqnarray}
or
\begin{eqnarray}
L(\nu)&=& \frac{8 \pi h \cos i}{c^2} \int_{r_i}^{r_f} \int_{0}^{2\pi}  \frac{\nu_e^3 r d \phi dr}{\exp{(\frac{h \nu_e}{k_{\rm B} T})}-1},
\end{eqnarray}
in this equation $i$ is the thin accretion disk's inclination angle around the black hole while $d$ is the distance between the observer and the centre of the thin accretion disc, along with the Planck constant $h$, emission frequency $\nu_e$, as well as the Boltzmann constant $k_{\rm B}$. Furthermore we need to the red-shift factor which can be found from
\begin{eqnarray}
g=\frac{\nu}{\nu_e}=\frac{k_\mu u^\mu_o}{k_\mu u^\mu_e},
\end{eqnarray}
where $\nu$ is the radiation frequency in the faraway observer's local rest frame. We also have the observer's four-velocity,
\begin{eqnarray}
u^\mu_o=\left(1,0,0,0\right),
\end{eqnarray}
and the the 4-velocity of the emitter
\begin{eqnarray}
u^\mu_e=\left(u^t_e,0,0,\Omega u^t_e\right).
\end{eqnarray}
We can take $r_i=r_{\rm ms}$ and $r_f$ to any arbitrary large distance from the compact object. To illustrate the effect, and see how the parameter $l$ affects the emission spectrum, we calculate the radiation spectrum $\nu L(\nu)$ numerically. We find that with the increase of the parameter $l$ for both accretion rates the observed luminosity increases as shown in Fig. 8. If we take the accretion case dot $M_0 \sim 8\,\times 10^{-5} M_{\odot} /{\rm yr}$, we find that the maximal value for the luminosity at particular $l$ occurs at the frequency $\nu=4.33 \times 10^{13}\, {\rm Hz}$ along with the following values for the luminosity:
\begin{equation}\notag
l=0.0:\,\nu L_{max}\simeq 1.158 \times 10^{33} {\rm erg}\; s^{-1},
\end{equation}
\begin{equation}\notag
l=1.8:\,\nu L_{max}\simeq 1.169 \times 10^{33} {\rm erg}\; s^{-1},
\end{equation}
\begin{equation}\notag
l=3.0:\,\nu L_{max}\simeq 1.188 \times 10^{33} {\rm erg}\; s^{-1}.
\end{equation}

Finally, we can think about accretion efficiency, which is calculated as the ratio of the rate of photons escaping from the disc surface to infinity compared to the rate of mass-energy transport to the black hole. There is a simple way to estimate the efficiency $\epsilon$ of the accretion disk, since it can be found to be proportional to the specific energy of the moving particle in the disc, as measured at the marginally stable orbit by
\begin{eqnarray}
\epsilon = 1- \tilde E_{\rm ms}.
\end{eqnarray}

\begin{figure}
    \centering
    \includegraphics[scale=0.63]{eps1.pdf}
    \caption{The accretion efficiency $\epsilon$ which is shown to be constant and does not depend on $l$ for $M=1$.}
\end{figure}

Interestingly, for the efficiency of the accretion disk in the Simpson-Visser spacetime we find that this quantity is not affected by the parameter $l$ and it is given by
\begin{eqnarray}
\epsilon = \frac{3-2 \sqrt{2}}{3}\simeq 0.0571909588.
\end{eqnarray}
Although the radiating energy changes, we found that based on the accretion efficiency $\epsilon$ one cannot distinguish the Simpson-Visser metric from the Schwarzschild geometry. This quantity is important and provides information about the efficient engine, in our case the accretion disk, for converting accreting matter's energy into electromagnetic radiation. Its interesting to see that this quantity remains constant and does not depend on the parameter $l$. This fact can be seen from Fig. 9.

\section{Conclusions}
The Simpson-Visser model has the additional benefit of analytically smoothly shifting from black holes (singular/regular) to wormholes (one-way/traversable). In some ways, the regular black hole shown above in this paper is in fact, the conventional Schwarzschild solution for $l = 0$. This is a clear contrast to the scenario in which a collapsing regular black hole “bounces” back into our own universe, and it is a situation worth considering in its own right.  It'll be interesting to see if we can tell the difference between a regular and singular black hole and a wormhole. This metric is particularly interesting because, depending on the value of the regularisation parameter $l$, it can either describe a regular/singular black hole or a traversable/one-way wormhole.
Therefore, in this paper, we studied the optical appearance of a thin accretion disc in the Simpson-Visser spacetime.  We calculate and illustrate the observed flux distributions as viewed by a distant observer. The image generated in the Simpson-Visser metric looks similar to that of Schwarzschild. In addition, we studied the electromagnetic properties of the accretion disk such as the temperature and the radiation flux of the energy by the accretion disk. It is shown that by increasing the regularized parameter $l$, the temperature, as well as the radiating flux increases, compared to the Schwarzschild geometry. Thus, for these solutions, we conclude that
the specific signatures, that appear in the electromagnetic spectrum, lead to the possibility of distinguishing
wormhole geometries from the Schwarzschild solution by
using astrophysical observations of the emission spectra
from accretion disks. As an interesting result, we show that the accretion efficiency remains constant and does not depend on $l$. In the near future, we plan to use astrophysical data to constrain the regularized parameter $l$.


\section{Introduction}
\label{sect:intro}
The attitude toward cryptocurrencies has radically changed in recent years: the initial skepticism that once accompanied them has given way to enthusiastic pronouncements, as their market capitalization reaches new record levels. Nonetheless, Bitcoin share of the entire cryptocurrency capitalization is above 40\% and concern is rising that it may become the ultimate speculative asset class, as people and financial institutions get increasingly involved with it.\footnote{According to \citet{coinmarketcap}, at the time of writing (September 2021) the global crypto market cap amounts to \$2.11 trillion, more than the current GDP of countries such as Italy and Canada \citep{worldbank}. For a detailed analysis on the growing interactions between Bitcoin and the financial sector, see \citet{bloombergcryptos}.}

A turning point in Bitcoin history could have been the launch of the first two regulated Bitcoin futures
by the Chigago Board of Exchange (CBOE) and the Chicago Mercantile Exchange (CME) in December 2017. Before that time, Bitcoin derivatives were only traded over-the-counter but since the introduction of a regulated derivative market, the interest of the financial sector toward Bitcoin has surged: hedge funds are investing in Bitcoin and banks are offering safekeeping services to their clients \citep{Leask:2020, Ossinger:2021, Sullivan:2021}.

Whether this is enhancing price discovery or increasing volatility of an already hyper-volatile asset is still unclear and the literature is far from a consensus. \citet{Hafner:2020} find evidence of a bubble-like behavior in Bitcoin returns and growing volatility in late 2017, suggesting it might be due to the new futures contracts and the increased media coverage during this period. Conversely, \citet{Shi:2017} finds a significant decrease in the spot volatility after the launch of CBOE contract, whereas later studies show that Bitcoin volatility increased right after the launch of CME future \citep{Kim:Lee:Kang:2020} as well as around its announcement date \citep{Corbet:Lucey:Peat:Vigne:2018}. Interestingly, the peak reached by Bitcoin price on December 16, 2017 matched the launch of the CME future; since then and for over a year, prices have experienced a sharp decline, suggesting that the newly introduced instruments may have allowed the entrance of ``pessimistic’’ traders willing to bet for a price drop but unable to do so until the creation of a derivative market \citep{Hale:Krishnamurthy:Kudlyak:Shultz:2018}.

Albeit the relatively few studies on the effect of Bitcoin futures, the impact of derivatives trading on the underlying spot volatility has been thoroughly investigated for instruments such as stocks and financial indexes. In particular, there are two conflicting theories about the effect of futures on the underlying spot markets and empirical evidence is mixed: some warn that speculative behaviors and information asymmetries brought by futures markets may increase volatility and jeopardize the value formation process in the spot market \citep{Stein:1987,Figlewski:1981,Antoniou:Holmes:1995,Harris:1989}; others argue that futures trading enhances information flows and price discovery, thereby reducing volatility and stabilizing the underlying spot market prices \citep{Danthine:1978,Moriarty:Tosini:1985,Edwards:1988,Bessembinder:Seguin:1992,Antoniou:Holmes:Priestley:1998,Jochum:Kondres:1998}. Recent works on individual stock futures \citep{Mckenzie:Brailsford:Faff:2001} and emerging economies \citep{Baklaci:Tutek:2006, Chen:Han:Li:2013,Bohl:Diesteldorf:Siklos:2015} support the idea that futures trading acts as a stabilizing force.

Contributing to this debate, in this paper we estimate the causal effect of both the CBOE and the CME instruments on Bitcoin volatility; we also investigate the presence of possible announcement effects and the causal impact on Bitcoin trading volume.

Indeed, the positive association between volume and volatility is well documented in the literature (e.g.,\citet{Bessembinder:Seguin:1993, Chen:Firth:Rui:2001, Kim:2005, Gebka:Wohar:2013}); therefore, part of the positive (negative) effect that we might observe on Bitcoin volatility could be the result of an indirect impact exerted by increased (decreased) trading volume. To verify this claim, we also estimate the effect of Bitcoin futures on trading volume and we explore the presence of a possible connection to volatility. Recent studies have attempted to investigate whether the positive relation between volume and volatility is present in the Bitcoin market as well, but the evidence is still mixed. Using squared returns as a proxy of volatility and a methodology based on non parametric causality-in-quantiles test, \citet{Balcilar:Bouri:Gupta:Roubaud:2017} find that trading volumes cannot predict the volatility of Bitcoin returns at any point of the conditional distribution. However, opposite results are obtained by the authors when repeating the analysis using GARCH volatility. Similarly, \citet{Bouri:Lau:Lucey:Roubaud:2019} show that trading volume do not Granger-cause Bitcoin volatility as proxied by squared returns. Conversely, using realized volatility computed from high-frequency data, \citet{Aalborg:Molnar:deVries:2019} find that volumes Granger cause volatility.

From the methodological viewpoint, most of the aforementioned studies employ Granger causality and only \citet{Kim:Lee:Kang:2020} uses a causal approach based on counterfactual outcomes --- Difference-in-Difference (DiD)--- in the attempt to attribute the uncovered effect to the introduction of futures.
Nonetheless, DiD relies on the assumption that the treated unit in the absence of intervention would have experienced the same trend as the control unit; therefore, this approach is impractical in all situations where, as in our case, there is no reliable control group available or the parallel trend assumption is troubled \citep{ONeill:Kreif:Grieve:Sutton:Sekhon:2016}.

In this paper, we employ the novel C-ARIMA approach proposed in \citet{Menchetti:Cipollini:Mealli:2021}, which, in contrast to DiD, allows to estimate causal effects in time series settings where suitable control series are unavailable and without imposing parallel trends. Furthermore, unlike traditional intervention analysis \citep{Box:Tiao:1975,Box:Tiao:1976, Bhattacharyya:Layton:1979, Larcker:Gordon:1980, Balke:Fomby:1994}, this method enables the computation of properly defined ``causal'' effects.
Indeed, C-ARIMA is based on the Rubin's potential outcomes approach \citep{Rubin:1974, Rubin:1975, Rubin:1978, Imbens:Rubin:2015} that allows to define the effect of an intervention as a contrast of ``potential outcomes'' as well as to discuss assumptions enabling the attribution of the effect to the intervention. Following a causal framework clarifying the assumptions needed to define and estimate the effect from available data, C-ARIMA uses ARIMA models to estimate the causal impact of the intervention by contrasting the observed post-intervention series with the counterfactual series that we would have had in the absence of intervention.

Although it was developed in a setting where a single intervention occurs, in this paper we extend the C-ARIMA approach to estimate causal effects in a multi-intervention setting, namely, a situation where a time series is subject to multiple interventions over time.
We then use it to evaluate the impact of the first two regulated Bitcoin futures on Bitcoin volatility as proxied by the Garman-Klass estimator and on transaction volumes. Our results indicate that CME future triggered an increase in both outcomes; there is also evidence that the positive association between volumes and volatility holds in the Bitcoin market as well and that the higher volumes induced by the CME instrument contributed, at least in part, to the observed effect on volatility. By controlling for these higher volumes, we find that Bitcoin volatility more than doubled due to the CME contract. Conversely, the CBOE instrument had only a minor impact on volumes and we do not find evidence of announcement effects.
Overall, our results seem to support the idea that futures trading is driving Bitcoin volatility upward, instead of diminishing it; this might be due to a renovated interest among investors following a perception of greater transparency, or to increased speculation, since regulated Bitcoin futures allow a greater number of people to bet against Bitcoin price. Regulators should be aware of these implications as banks and financial institutions are opening up to this new investment opportunity.

The rest of the paper is structured as follows: Section \ref{sect:background} describes the background for our empirical analysis; Section \ref{sect:causal_framework} sets out the causal framework; Section \ref{sect:cARIMA} extends the C-ARIMA model to the multi-intervention setting; Sections \ref{sect:data_method} and \ref{sect:results} illustrate the data and the results of the empirical analysis; Section \ref{sect:conclusions} concludes.

\section{Background}
\label{sect:background}

Bitcoin is a peer-to-peer payment system created in 2009
under the pseudonym of Satoshi Nakamoto \citep{Nakamoto:2008}. In contrast to fiat money relying on central banks and intermediaries, Bitcoin is decentralized, meaning that its value is not backed by any central bank.
In short, transaction data are recorded in \textit{blocks}, each containing a reference (the hash) to the previous ones, thereby forming a chain know as the \textit{blockchain}. Transactions are validated through the so-called \textit{mining} process, which requires to find the hash by solving a time consuming cryptographic problem. Every time this happens, the new validated block is added to the chain and miners are rewarded with new Bitcoins. Thus, essentially, the blockchain constitutes an electronic public ledger of validated transactions, which is stored and updated on miner's computers (the network nodes) to avoid double spending problems and frauds. Mining is the only way new Bitcoins are introduced in the system and their supply is limited by design to a maximum of 21 million units that will be reached by 2140.\footnote{More precisely, an algorithm adjusts the difficulty of the numerical problem for finding the hash based on the network performance, so that new blocks are generated at a fixed rate (every 10 minutes); in addition, the number of generated Bitcoins is halved by design every 210,000 blocks (approximately 4 years). Considering these facts, it can be estimated that by 2140 miners' reward will be roughly zero, meaning that the maximum supply will be reached.}

Bitcoins are traded on multiple exchanges, the top A-rated ones being Coinbase, Gemini and Bitstamp \citep{CryptoBenchmark:2021}\footnote{In the Exchange Benchmark Report released in February 2021, CryptoCompare rated 162 exchanges according to 68 qualitative and quantitative metrics. They show that existing metrics, such as volume or liquidity, can be easily manipulated (e.g. volumes can be inflated through strategies such as trading competitions, airdrops and transaction fee mining) and thus are inadequate to reflect the reliability of the trading venue. For further details regarding these strategies and the rating methodology see \citet{CryptoBenchmark:2021}.}. At the time of writing (September 2021, timestamp 2021-09-16 00:00:00 GMT) Bitcoin is trading at \$47,629.49 while in 2010 1BTC was valued \$0.05.
Such extreme price fluctuations (and the related opportunities of huge profits) probably contributed to spark investor interest toward cryptocurrencies and
to the decision, by two major derivative exchanges, of starting a regulated derivative market for Bitcoin futures. Specifically, on August 2, 2017 CBOE announced its partnership with
Gemini Trust Company to use Gemini’s Bitcoin market data in the creation of derivatives products for trading; the future was then released on December 10. Meanwhile, CME  announced the launch of its contract becoming effective starting from December 18.
Since then, having seen its market share quickly eroded by the new future \citep{Baydakova:2019}, CBOE has stopped listing additional Bitcoin futures for trading and after nine months from the release of the first contract, also the Intercontinental Exchange (ICE) started offering Bitcoin derivatives. Table \ref{tab:event_history} summarizes the major events occurring in Bitcoin futures' history up to December 2017.

In our empirical analysis, we focus on the futures launched by CBOE and CME, as it is commonly done in the literature, to assess whether regulated futures have a role in determining Bitcoin volatility and volume.

To improve the prediction accuracy of our models, we include covariates that might be linked to Bitcoin trading.
Early works evidence association between Bitcoin price and variables such as the exchange-trade ratio (i.e. ratio between trade volume and transaction volume), hash rate, Google and Wikipedia queries and the Shangai index \citep{Kristoufek:2013,Kristoufek:2015,Bouoiyour:Selmi:2015}. In a later study,  \citet{Xin:Wang:2017} find a strong association between Bitcoin price and key economic fundamentals (US interest rate, USD money supply, Bitcoin money supply and transaction volume); they also find a connection between technology factors such as mining difficulty and public recognition and the price before 2014.\footnote{A major event that affected the Bitcoin network is the failure of Mt. Gox, the leading Bitcoin exchange until February 28, 2014, when it filed for bankruptcy protection after announcing a theft of about 850,000 BTC following a security breach \citep{McLannahan:2014,Cermak:2017}. This is considered in the literature as a path-breaking event and many works focusing on Bitcoin returns make a distinction between an early period (before Mt. Gox failure) and a later period (after the failure).} Similarly, \citet{Ciaian:Kancs:Rajcaniova:2016} evidence that market forces of Bitcoin supply and demand are strongly related to Bitcoin price changes while proxies of public recognition (i.e. views on Wikipedia, new posts and new members on Bitcoin forums) are mainly associated with the price before 2014. Finally, \citet{Liu:Tsyvinski:2021} find that network factors such as the number of wallet users and the number of active addresses, are important drivers of cryptocurrency price. In Section \ref{sect:data_method} we describe how to include these drivers in our models.

\begin{table}[!h]
\centering
\caption{Major events in Bitcoin futures' history up to December 2017}
\label{tab:event_history}
\begin{tabular}{ll} \toprule
Date & Event \\ \midrule
2017/08/02 & CBOE announces launch of futures by 2017-Q4 \\
2017/10/31 & CME announces launch of futures by 2017-Q4\\
2017/12/01 & CME announces launch of futures on Dec 18 \\
2017/12/04 & CBOE announces launch of futures on Dec 10 \\
2017/12/10 & Launch of CBOE futures \\
2017/12/18 & Launch of CME futures \\
\bottomrule
\end{tabular}
\end{table}

\section{Causal Framework}
\label{sect:causal_framework}

In this section we illustrate the framework that enables us to make causal statements in a multi-intervention time series setting where the treatments are not randomized.  We first discuss the assumptions on the treatments, potential outcomes and covariates, then we define the causal estimands. In doing so, we refer to the causal framework that is commonly adopted in observational settings with a single intervention \citep{ Choirat:Dominici:Mealli:Papadogeorgou:Wasfy:Zigler:2018, Menchetti:Bojinov:2021,Menchetti:Cipollini:Mealli:2021} and, when needed, we integrate it in order to encompass the multi-intervention situation.

\subsection{Assumptions}
\label{subsect:assumptions}

For a generic statistical unit, let $\W_t \in \{ 0,1 \}$ be the treatment assignment at time $t \in \{1,\dots, T\}$, where $\{1\}$ denotes that a ``treatment'' (or ``intervention'') has taken place and $\{0\}$ denotes control. Then, $\W_{1:T} = (\W_1,\dots,\W_T) $ is the assignment path, i.e., the sequence of treatments received by the unit. A realization of $\W_t$ is denoted with the lower case letter $\w_t$.

In our empirical application, the statistical unit is Bitcoin cryptocurrency and the main goal of the analysis is investigating the effect on both volatility and volume generated by the launch of Bitcoin futures by two major regulated exchanges. Since both exchanges disclosed their plans to develop Bitcoin futures, we have two types of interventions: i) the announcements about the upcoming futures; ii) the actual introduction of the two contracts. To make reliable estimates of the effects attributable to these interventions, in accordance with the causal framework presented in \citet{Menchetti:Cipollini:Mealli:2021} we need to discuss three assumptions: non-anticipating potential outcomes, covariates-treatment independence and non-anticipating treatments.

When a statistical unit is assigned to treatment, its outcome may differ from the one we would have observed if the same unit had been assigned to control. In the Rubin causal framework these are called ``potential outcomes''. Let $\Y_t (\w_t)$ denote the potential outcome at time $t$ for a treatment $\w_t$: \textit{non-anticipating potential outcomes} can be function of past treatments but they are unaffected by future treatments and the dependency on the past treatment sequence is denoted by $\Y_t(\w_{1:t})$. In our application, we admit that today volatility is the result of previous announcements and of the trading activity of existing futures but we exclude any influence arising from future announcements and future contracts. In other words, we are ruling out the possibility that market participants have access to privileged information.

In case of binary treatments, like the ones we are considering, at each point in time there are two potential outcomes corresponding, respectively, to the outcome in the presence and in the absence of the intervention, i.e., $\Y_t(1)$ and $\Y_t(0)$. When the intervention occurs, its causal effect is given by direct comparison between them; however, only $\Y_t(1)$ is actually observed, whereas $\Y_t(0)$ is commonly referred as the ``missing'' or ``counterfactual'' outcome and needs to be estimated. In this process, informative covariates may help.

As outlined in Section \ref{sect:background}, Bitcoin price might be related to several economic and technology factors, such as USD money supply and mining difficulty. Including their values in our analysis may enhance the prediction accuracy of the counterfactual outcome and, as a result, the reliability of the estimated causal effect. However, we should exclude from our covariates set any predictor that might be influenced by the intervention, otherwise the estimated effect will be biased. This \textit{covariates-treatment independence} assumption also applies in our multi-intervention setting; therefore, in the empirical application we use as predictors only those covariates for which such assumption is plausible. For example, since the mining process is ruled by an algorithm and new Bitcoins are created at a fixed rate, the current Bitcoin supply is unaffected by the introduction of futures; conversely, network factors such as the number of wallets and all the proxies of investors' interests (e.g., Google and Wikipedia queries) may have increased following the introduction of the new futures and thus can not be assumed independent of the treatment.

We also assume \textit{non-anticipating treatments}, meaning that the assignment mechanism (i.e., the process that determines the units receiving treatment) depends solely on past outcomes and past covariates. To understand the importance of such an assumption, consider as an example a research on a new drug where doctors are asked to select the participants among their patients. If doctors assign to treatment only those patients they believe have better chance to complete the treatment successfully and without side effects, the results of such study would be biased due to a ``confounded'' assignment to treatment. A non-anticipating treatment in a time series setting is the analogous of the unconfounded assignment mechanism in the cross-sectional setting \citep{Imbens:Rubin:2015}, since it
is essential in ensuring that, conditioning on past outcomes and covariates, any difference in the potential outcomes is due to the treatment. As a result, in our application, the effect of the CBOE future is correctly defined up to one week after its launch (beyond that point, the effect would be confounded by the introduction of the CME contract); meanwhile, this assumption ensures that the effect of the CME contract is not confounded by the first future, since we condition on past events.

In a time series setting, interventions can occur at every time point, like the announcements made by the two exchanges. In this case, the total number of potential outcomes at the end of the analysis period is $2(2^T-1)$, corresponding to $2^T$ potential paths (Figure \ref{fig:potential_paths} provides an illustration in a simple multi-intervention setting). However, it is not uncommon to observe a single intervention affecting the time series in a persistent way, such as a price discount for several weeks in a row during a marketing campaign \citep{Brodersen:Gallusser:Koehler:Remy:Scott:2015, Menchetti:Cipollini:Mealli:2021, Menchetti:Bojinov:2021}. As in the case of the two futures, we may also have multiple persistent interventions during the analysis period.

\begin{assumpt}[\textit{N} Persistent interventions]
\label{assumpt:n_persistent_int}
Indicating with $\Lambda = \{t_1,\dots,t_N \}$ the subset of time points at which the interventions take place, we say that the unit received $N$ persistent interventions, if for all $t \le t_1$ we have $\W_{t}=(0,\dots,0)$ and for all $t,t' \in \{ t_{n},\dots,t_{n+1}-1\}$ we have $\W_{t}=\W_{t'}$; we also assume that the effects generated by the N persistent interventions are additive.
\end{assumpt}

In words, at each time point in the interval between two subsequent interventions, the statistical unit is either always treated or always assigned to control. Additivity of the effects follows quite naturally: if the interventions persist over time, their effect will likely do the same. Therefore, the observed outcome is the result of all previous interventions.\footnote{Notice that the assumption of additive effects does not pose any strong restriction: in most cases, we can still recover an additive structure by applying simple transformations (e.g., logarithmic transformation).}
In our empirical application, the two futures qualify as persistent interventions, since they are bound to affect permanently the dynamics of Bitcoin price. Indeed, even though the exchange can withdraw its future from the market at any time, the future would trade until its expiration date, which is standardized and set up in advance.
In this special setting, Assumption \ref{assumpt:n_persistent_int} narrows down the set of potential paths: instead of $2^T$ potential paths, in the time interval between two persistent interventions, $ t \in \{t_{n},\dots, t_{n+1}-1 \}$ we only have two possible paths,
$$\w_{t_{n}:t_{n+1}-1} = (\underbrace{1,\dots,1}_{t \in \{t_n,\dots,t_{n+1}-1\}})\hspace{8pt}
\text{and} \hspace{8pt} \w'_{t_{n}:t_{n+1}-1} = (\underbrace{0,\dots,0}_{t \in \{t_n,\dots,t_{n+1}-1\}})$$
where $\w_{t_{n}:t_{n+1}-1} = \bs{1}$ is a vector of the observed treatment path indicating that the $n$-th persistent intervention has occurred and persists over time, $\w'_{t_{n}:t_{n+1}-1} = \bs{0}$ is the counterfactual treatment path. Figure \ref{fig:potential_paths_persistent} gives an example of the potential outcomes and treatment paths in a simple case of a single persistent intervention.

\begin{figure}[!h]
\centering
\caption{Potential outcome time series when $T=3$ and two treatments are administered at times $t_1 = 1$ and $t_2=3$. The total number of potential outcomes is $2(2^3-1) = 14$ corresponding to $2^3$ paths. Only one path is actually observed (indicated by the solid line) whereas the others are missing (indicated by the dashed lines).}
\label{fig:potential_paths}
\begin{tikzpicture}
\tikzstyle{level 1}=[level distance=3cm, sibling distance=3cm]
\tikzstyle{level 2}=[level distance=3.5cm, sibling distance=1.5cm]
\tikzstyle{level 3}=[level distance=3.5cm, sibling distance=0.6cm]
\node{Y}[grow=right]
 child[dashed]{node{$\Y_1(0)$}
  child{node{$\Y_2(0,0)$}
   child{node{$\Y_3(0,0,0)$}}
   child{node{$\Y_3(0,0,1)$}}
  }
  child{node{$\Y_2(0,1)$}
   child{node{$\Y_3(0,1,0)$}}
   child{node{$\Y_3(0,1,1)$}}
  }
 }
 child{node{$\Y_1(1)$} edge from parent[black]
  child[black]{node{$\Y_2(1,0)$} edge from parent[black]
   child[dashed]{node{$\Y_3(1,0,0)$}}
   child{node{$\Y_3(1,0,1)$}}
  }
  child[dashed]{node{$\Y_2(1,1)$}
   child[black]{node{$\Y_3(1,1,0)$}edge from parent[black]}
   child{node{$\Y_3(1,1,1)$}}
  }
 };
\end{tikzpicture}
\end{figure}

\begin{figure}[h!]
\centering
\caption{Potential outcome time series of a single persistent intervention occurring at time $t = 3$. The solid line represents the observed series while the dashed line depicts the missing potential outcome series.}
\label{fig:potential_paths_persistent}
\begin{tikzpicture}
\tikzstyle{level 1}=[level distance=3cm, sibling distance=3cm]
\tikzstyle{level 2}=[level distance=3.5cm, sibling distance=1.5cm]
\tikzstyle{level 3}=[level distance=3.5cm, sibling distance=1.5cm]
\node{Y}[grow=right]
 child{node{$\Y_1(0)$}
  child{node{$\Y_2(0,0)$}
   child[dashed]{node{$\Y_3(0,0,0)$}
    child[dashed]{node{$\Y_4(0,0,0,0)$}}}
   child{node{$\Y_3(0,0,1)$}
    child{node{$\Y_4(0,0,1,1)$}}}
  }
 }
 ;
\end{tikzpicture}
\end{figure}

\subsection{Causal Estimands}
\label{subsect:estimands}

We can now define the causal estimands of interest. In particular, we provide definitions for two causal effects: the contemporaneous and the pointwise causal effects. The former is used to describe the effect of the announcements, which is reasonably short-term, whereas the latter is used to describe the impact produced by the futures, which qualify as persistent interventions.

\begin{defn}[Contemporaneous effects]
Indicating with $\Lambda = \{t_1, \dots, t_N \}$ the subset of time points at which the active treatment is administered, the contemporaneous causal effect of the $n$-th treatment at time $t_n \in \Lambda$ conditioning on the observed treatment path $\w_{1:t_n-1}^{obs}$ is,
\begin{equation}
\label{eqn:cont}
\tau^{(n)}(\w_{1:t_n-1}^{obs}, 1;0) = \Y_{t_n}(\w_{1:t_n-1}^{obs},1) - \Y_{t_n}(\w_{1:t_n-1}^{obs},0) .
\end{equation}
\end{defn}
In words, this is the instant effect of a treatment conditioning on the observed treatment path up to that time point. The contemporaneous effect is similar to the general effect in \citep{Bojinov:Shephard:2019}, with the difference that the latter estimand refers to an experimental setting where the intervention is randomly allocated at any point in time and thus it is not necessary to condition on the treatment path.

\begin{example}
Assume that the solid line in Figure \ref{fig:potential_paths} represents the observed path of the announcements. At $t_1=1$ we immediately observe an intervention and, since there is no past history, the contemporaneous effect is $\tau^{(1)}(1;0) = \Y_1(1) - \Y_1(0)$. At $t_1+1=2$ we observe control. Then, at $t_2=3$ we have a second announcement and by conditioning on the observed path, we can define its contemporaneous causal effect as $\tau^{(2)}(\w^{obs}_{1:2}, 1;0) = \Y_3(1,0,1) - \Y_3(1,0,0)$.
\end{example}

We now define a class of causal estimands measuring the effect of the $n$-th persistent intervention. Recall that, conditioning on previous treatments, there are only two possible potential paths in the time interval between two persistent interventions, $\w_{t_n:t_{n+1}-1} = \bs{1}$ and $\w'_{t_n:t_{n+1}-1} = \bs{0}$.

\begin{defn}[Pointwise effects]
Indicating with $\Lambda = \{t_1, \dots, t_N \}$ the subset of time points at which $N$ persistent interventions take place, the point effect at time $t \in \{t_n,\dots,t_{n+1}-1 \}$ of the $n$-th persistent intervention conditioning on the observed path $\w^{obs}_{1:t_n-1}$ is
\begin{equation}
\label{eqn:point}
\tau_t^{(n)}(\w^{obs}_{1:t_n-1},\bs{1}; \bs{0} ) = \Y_t(\w^{obs}_{1:t_n-1}, \bs{1}) - \Y_t(\w^{obs}_{1:t_n-1}, \bs{0}).
\end{equation}
Thus, the cumulative pointwise effect of the $n$-th intervention up to time $t$ is,
\begin{equation}
\label{eqn:cum_pointwise}
\Delta_t^{(n)}(\w^{obs}_{1:t_n-1},\bs{1}; \bs{0}) = \sum\limits_{s=t_n}^t\tau_s^{(n)}(\w^{obs}_{1:t_n-1},\bs{1}; \bs{0})
\end{equation}
and the temporal average pointwise effect of the $n$-th intervention is,
\begin{equation}
\label{eqn:temp_avg_pointwise}
\bar{\tau}_t^{(n)}(\w^{obs}_{1:t_n-1},\bs{1}; \bs{0}) = \frac{\Delta_t^{(n)}( \w^{obs}_{1:t_n-1},\bs{1}; \bs{0})}{t-t_n + 1} .
\end{equation}
\end{defn}
Notice that for $t = t_n$, the pointwise effect collapses to the contemporaneous effect and in case of one single persistent intervention ($N = 1$) it matches the pointwise effect as defined in \citet{Menchetti:Cipollini:Mealli:2021}.

\begin{example}
Assume that CBOE introduces Bitcoin futures at time $t_1=3$, as outlined in Figure \ref{fig:potential_paths_persistent}. The point effect at $t_1$ is the contemporaneous effect $\tau^{(1)}(\w^{obs}_{1:2},1;0) = \Y_3(0,0,1) - \Y_3(0,0,0)$ and the point effect at time $t_1+1 = 4$ is $\tau_{4}^{(1)}(\w^{obs}_{1:2},1,1;0,0) = \Y_4(0,0,1,1) - \Y_4(0,0,0,0)$. Then, the cumulative and the temporal average pointwise effects are, respectively, $\Delta_4^{(1)}(\w^{obs}_{1:2},1,1;0,0) = \tau^{(1)}(\w^{obs}_{1:2},1;0) +  \tau_{4}^{(1)}(\w^{obs}_{1:2},1,1;0,0)$ and $\bar{\tau}_4^{(1)}(\w^{obs}_{1:2},1,1;0,0) = \frac{1}{2}\Delta_4^{(1)}(\w^{obs}_{1:2},1,1;0,0)$.
\end{example}

Thus, the contemporaneous effect can also be interpreted as a limiting case of the pointwise effect, namely, the effect of a ``pulse'' intervention that lasts only for one time point. As a result, in the next section, the estimators of the contemporaneous effects will be derived as limiting cases of the pointwise effects.

\section{C-ARIMA}
\label{sect:cARIMA}
In this section, we extend the C-ARIMA model \citep{Menchetti:Cipollini:Mealli:2021} to the multi-intervention setting and we use it to derive estimators for the causal effects defined in Section \ref{subsect:estimands}.

\subsection{The model}

The general formulation of the C-ARIMA model in a setting with a single intervention occurring at time $t^*$ is a linear regression with seasonal ARIMA errors and the addition of a component $\tau_t$,
\begin{equation}
\label{eqn:arima_general}
(1-L^s)^D (1-L)^d \Y_t(\w)  =  \frac{\Theta_Q (L^s)\theta_q(L) }{ \Phi_P(L^s) \phi_p(L)}\varepsilon_t +  (1-L^s)^D (1-L)^d \X_t' \beta + \tau_t \mathds{1}_{\{\w = 1\}}
\end{equation}
where $(1-L^s)^D$ and $(1-L)^d$ are the differencing operators; $\theta_q(L)$ and $\phi_p(L)$ are lag polynomials having roots all outside the unit circle; $\varepsilon_t$ is white noise with mean 0 and variance $\sigma^2_{\varepsilon}$; $\X_t$ is a set of covariates satisfying the covariates-treatment independence assumption; $\Theta_Q (L^s)$, $\Phi_P(L^s)$ are the lag polynomials of the seasonal part of the model with period $s$ and having roots all outside the unit circle. More importantly, $\tau_t = 0$ $\forall t < t^*$ and $\mathds{1}_{\{\w = 1\}}$ is an indicator function which is 1 if the statistical unit receives the treatment. As highlighted in \citet{Menchetti:Cipollini:Mealli:2021}, this means that $\tau_t$ can be interpreted as the causal effect of the treatment at time point $t$, since it is defined as a contrast of potential outcomes. To see that, indicate with $T(\cdot)$ the transformation of $\Y_t(\w)$ needed to achieve stationarity, i.e. $T(\Y_t(\w)) =(1-L^s)^D (1-L)^d \Y_t(\w) $. Notice that the same transformation is also applied to $\X_t$. Then, we can define
\begin{equation*}
z_t = \frac{\Theta_Q (L^s)\theta_q(L) }{ (1-L^s)^D (1-L)^d \Phi_P(L^s) \phi_p(L)}\varepsilon_t
\end{equation*}
so that model (\ref{eqn:arima_general}) becomes,
\begin{equation*}
\Y_t(\w) =  z_{t} + \X_{t}' \beta + \tau^Y_{t} \mathds{1}_{\{\w = 1\}}
\end{equation*}
where
\begin{equation*}
\tau^Y_t = \frac{\tau_t}{(1-L^s)^D (1-L)^d}
\end{equation*}
is the causal effect on the untrasformed variable, defined as a contrast of potential outcomes,
i.e., $\tau^Y_t = \Y_t(\w = 1) - \Y_t(\w = 0)$. Similarly, in a multi-intervention setting where $N$ persistent interventions take place at time points $\Lambda \in \{ t_1,\dots,t_N\}$, under non-anticipating treatments and non-anticipating potential outcomes, the outcome at time $t \in \{t_n,\dots,t_{n+1}-1 \}$ is the result of all past interventions,
\begin{equation}
\label{eqn:arima_general_multi}
\Y_{t}(\w) = z_{t} + \X_{t}' \beta + \sum\limits_{j = 1}^{n-1} \tau_{t}^{(j)} + \tau_{t}^{(n)} \mathds{1}_{ \{ \w = 1\} }
\end{equation}
where the summation comes from the additivity of the effects produced by $n-1$ persistent interventions and $\tau_t^{(n)}$ is the point effect at time $t$ of the $n$-th intervention.  Now, assume we observe $\Y_t$ up to time $t_n-1$ and that $\w = 0$. Then, for a positive integer $k$, the $k$-step ahead forecast of $\Y_t$ in the absence of the $n$-th intervention, given all the information up to time $t_n -1$ is,
\begin{equation}
\label{eqn:1step_forecast}
\hat{\Y}_{t_n-1+k}(0) = E[\Y_{t_n-1+k}(0) | \info{I}{t_n-1}]  = \hat{z}_{t_n-1+k|t_n-1}  + \X_{t_n-1+k}' \beta + \sum\limits_{j = 1}^{n-1} \tau_{t_n-1+k}^{(j)}.
\end{equation}
Notice that $\hat{\Y}_{t_n -1 +k}(0)$ is, by definition, the expectation of the outcome series when the $n$-th intervention does not occur;
thus, it can be considered an estimate of the missing potential outcomes at time $t_n -1 +k$ for a persistent intervention occurring at time $t_n$. Therefore, the point effect of the $n$-th intervention can be estimated by the quantity
\begin{align}
\label{eqn:pred_error}
\nonumber
\hat{\tau}^{(n)}_{t_n-1+k} & = \Y_{t_n -1 + k}(1) - \hat{\Y}_{t_n-1+k}(0) \\ \nonumber
& = z_{t_n -1 + k} + \sum\limits_{j=1}^{n-1} \tau_{t_n -1+k}^{(j)} + \tau_{t_n -1+k}^{(n)} - \hat{z}_{t_n - 1 + k|t_n-1} - \sum\limits_{j = 1}^{n-1} \tau_{t_n -1+k}^{(j)} \\
& = z_{t_n -1 + k} - \hat{z}_{t_n - 1 + k|t_n-1} + \tau_{t_n -1+k}^{(n)}.
\end{align}
Furthermore, the limiting case $k = 1$ leads to the $1$-step ahead forecast $\hat{\Y}_{t_n}$, which is an estimate of the missing potential outcome in the absence of intervention for the contemporaneous effect.

\subsection{Causal effect inference}
\label{subsect:estimators}
Notice that, since Equations (\ref{eqn:1step_forecast})-(\ref{eqn:pred_error}) are defined by conditioning on the information set up to time $t_n-1$, more formally we have $\hat{\Y}_{t_n -1 +k}(0) \equiv \hat{\Y}_{t_n -1 +k}(\w_{1:t_n-1}^{obs}, \bs{0})$ and $\Y_{t_n -1 +k}(1) \equiv \Y_{t_n -1 +k}(\w_{1:t_n-1}^{obs}, \bs{1})$. Based on the C-ARIMA model for the multi-intervention setting, we can now derive estimators for the pointwise and the contemporaneous causal effects.

\begin{defn}
Denote with $\Lambda = \{t_1,\dots,t_N \}$ the subset of time points at which $N$ persistent interventions take place and indicate with $\w_{1:t_n-1}^{obs}$ the observed treatment path. For a positive integer $k$, let $\Y_{t_n-1+k}(\w_{1:t_n-1}^{obs}, \bs{1})$ be the observed time series and let $\hat{\Y}_{t_n-1+k}(\w_{1:t_n-1}^{obs}, \bs{0})$ be the $k$-step ahead forecast as defined in (\ref{eqn:1step_forecast}). For any $k>0$, an estimator of the point  effect of the $n$-th intervention at time $t_n-1+k$ is,
\begin{equation}
\label{eqn:point_estimator}
\hat{\tau}_{t_n-1+k}^{(n)}(\w_{1:t_n-1}^{obs}, \bs{1}; \bs{0}) = \Y_{t_n-1+k}(\w_{1:t_n-1}^{obs}, \bs{1}) - \hat{\Y}_{t_n-1+k}(\w_{1:t_n-1}^{obs}, \bs{0}) .
\end{equation}
Then, an estimator of the cumulative pointwise effect of the $n$-th intervention up to time $t_n-1+k$ is,
\begin{equation}
\label{eqn:cumpoint_estimator}
\hat{\Delta}_{t_n-1+k}^{(n)}(\w_{1:t_n-1}^{obs},\bs{1}; \bs{0}) = \sum\limits_{h = 1}^{k} \hat{\tau}_{t_n-1+k}^{(n)}(\w_{1:t_n-1}^{obs},\bs{1}; \bs{0})
\end{equation}
and, finally, an estimator of the temporal average pointwise effect is,
\begin{equation}
\label{eqn:avgpoint_estimator}
\hat{\bar{\tau}}_{t_n-1+k}^{(n)}(\w_{1:t_n-1}^{obs},\bs{1}; \bs{0}) = \frac{\hat{\Delta}_{t_n-1+k}^{(n)}(\w_{1:t_n-1}^{obs},\bs{1}; \bs{0})}{k}.
\end{equation}
\end{defn}
Furthermore, by setting $k = 1$ in Equation (\ref{eqn:point_estimator}) we get an estimator of the contemporaneous effect of the $n$-th intervention at time $t_n$,
\begin{equation*}
\hat{\tau}^{(n)}(\w_{1:t_n-1}^{obs},1; 0) = \Y_{t_n}(1) - \hat{\Y}_{t_n}(0) .
\end{equation*}
Inference on the pointwise effects and, by extension, on the contemporaneous effect can be performed using hypothesis tests based on the above estimators. The following theorem illustrates their distributional
properties.

\begin{theorem}
Let $\{ \Y_t(\w)\}$ follow the regression model with ARIMA errors described by Equation (\ref{eqn:arima_general_multi}) and, for any $k >0$, let $H_0: \tau^{(n)}_{t_n-1+k} = 0$ the null hypothesis that the $n$-th intervention has no effect. Then, for some $b_i$ coefficients with $b_0 = 1$, the estimators of the point, cumulative and temporal average effects under $H_0$ can be expressed as,
\begin{align}
\label{eqn:point_estimator_var_h0}
\hat{\tau}_{t_n-1+k}^{(n)}(\w_{1:t_n-1+k}^{obs};\bs{1},\bs{0})|H_0 & = \sum\limits_{j = 1}^{k} \varepsilon_{t_n-1+j} \sum\limits_{i = 0}^{k-j} b_i \psi_{k-j-i} \\
\label{eqn:cumpoint_estimator_var_h0}
\hat{\Delta}_{t_n-1+k}^{(n)}(\w_{1:t_n-1+k}^{obs}, \bs{1}, \bs{0})|H_0 & = \sum\limits_{h = 1}^k \varepsilon_{t_n-1+h} \sum\limits_{i = 0}^{k-h} b_i \sum\limits_{j = i}^{k-h} \psi_{k-h-j} \\
\label{eqn:avgpoint_estimator_var_h0}
\hat{\bar{\tau}}_{t_n-1+k}^{(n)}(\w_{1:t_n-1+k}^{obs}, \bs{1}, \bs{0})|H_0 & = \frac{1}{k}\sum\limits_{h = 1}^k \varepsilon_{t_n-1+h} \sum\limits_{i = 0}^{k-h} b_i \sum\limits_{j = i}^{k-h} \psi_{k-h-j},
\end{align}
where the $\psi$'s are the coefficients of a moving average of order $k-1$ whose values are function of the ARMA parameters in Equation (\ref{eqn:arima_general_multi}). In case the error term $\varepsilon_t$ is assumed to be Normally distributed, Equations (\ref{eqn:point_estimator_var_h0})--(\ref{eqn:avgpoint_estimator_var_h0}) become,

\begin{align}
\label{eqn:point_estimator_var}
\hat{\tau}_{t_n-1+k}^{(n)}(\w_{1:t_n-1}^{obs},\bs{1};\bs{0})|H_0
& \sim
N \left[ 0, \sigma^2_{\varepsilon_{n-1}}\sum\limits_{j = 1}^{k} \left( \sum\limits_{i = 0}^{k-j} b_i \psi_{k-j-i} \right)^2 \right] \\
\label{eqn:cumpoint_estimator_var}
\hat{\Delta}_{t_n-1+k}^{(n)}(\w_{1:t_n-1+k}^{obs}, \bs{1}, \bs{0})|H_0
& \sim N \left[ 0,  \sigma^2_{\varepsilon_{n-1}}\sum\limits_{h=1}^k  \left(\sum\limits_{i = 0}^{k-h} b_i \sum_{j=i}^{k-h} \psi_{k-h-j} \right)^2 \right] \\
\label{eqn:avgpoint_estimator_var}
\hat{\bar{\tau}}_{t_n-1+k}^{(n)}(\w_{1:t_n-1+k}^{obs}, \bs{1}, \bs{0})|H_0
& \sim N \left[ 0, \frac{\sigma^2_{\varepsilon_{n-1}}}{k^2}\sum\limits_{h=1}^k  \left(\sum\limits_{i = 0}^{k-h} b_i \sum_{j=i}^{k-h} \psi_{k-h-j} \right)^2 \right],
\end{align}

where $\sigma_{\varepsilon_{n-1}}$ is the variance of the error term of the model estimated on the observations up to time $t_n-1$.

Proof: the derivation is analogous to Theorem 2 in \citet{Menchetti:Cipollini:Mealli:2021}, where we can replace the single intervention date $t^*$ with $t_n-1$, the day before the $n$-th intervention.

\end{theorem}

Thus, the estimation of the pointwise effect is performed in two steps.  First, we estimate a C-ARIMA model up to the day preceding the persistent intervention, i.e., the launch of the two futures by CBOE and CME. Then, we use the covariates and the estimated coefficients to forecast the outcome time series (i.e., Garman-Klass volatility proxy and trading volume) up to a pre-specified time horizon; for example, we may be interested in estimating the effect on volatility after few days or few months from the launch of Bitcoin futures. The difference between the observed and the predicted outcome is the estimated pointwise effect. Similarly, we can estimate the contemporaneous effects of each announcement by conditioning to the information set up to the day before the announcement and forecasting the Garman-Klass proxy and the trading volume one-step ahead. If the announcements or the launch of futures had an impact on Bitcoin volatility and volume, we would find a significant deviation from the forecasted outcomes.

\section{Empirical application}
\label{sect:data_method}
\subsection{Data}
Economic data for this analysis have been gathered from Bloomberg, while Bitcoin daily prices have been collected from CryptoCompare. Since Bitcoin is a (crypto)currency, its price is recorded in terms of another currency, USD in our case, meaning that the Bitcoin price is the BTC-USD exchange rate (number of USD needed to buy 1 BTC). The other Bitcoin related data (i.e. the hash rate and total number of Bitcoins in circulation) have been gathered from Blockchain.com.

Since Bitcoins are traded on multiple exchanges, the quotation of the BTC-USD rate is not unique, i.e. different exchanges trade Bitcoins at different prices. As explained in \citet{Cermak:2017}, the price varies across exchanges mainly because of different fee policies and cashout methods but those divergent standards and the slow verification process make arbitrage opportunities difficult to exploit\footnote{For example, on November 20, 2019 at 4.50 p.m the price for 1 BTC was \$8,104.39 on Coinbase and \$8,146.90 on Bitfinex, with a difference of \$42.51 (Source: CryptoCompare.com).}. Aside from the obvious economic implications, this means that the source we take our data from matters: exchange-based data providers present their own quotes, whereas external data providers (e.g. Bloomberg, CryptoCompare) compute their own index, usually a weighted average of all prices across major exchanges.

The goal of our analysis is estimating the effect of futures trading on Bitcoin volatility and volumes, and, since the upcoming futures were announced by several press releases, we also investigate possible announcement effects. We focus on the first two futures introduced by CBOE and CME, hence, the analysis period spans from May 2014 to January 2018. We start from May 2014 to avoid the market turbulence following the failure of Mt. Gox; ending the analysis period in January 2018 allows us to follow the Bitcoin network for three weeks after the launch of the CME future.

The BTC-USD daily volatility is proxied by the unbiased Garman-Klass estimator  \citep{Garman:Klass:1980,Molnar:2012}, computed as
\begin{align*}
\hat{\sigma}_{GK}^2 & = 0.5 (\ln(H)-\ln(L))^2 - (2\ln2 -1)(\ln(C) - \ln(O))^2 \\
\hat{\sigma}_{GK} & = \sqrt{\hat{\sigma}_{GK}^2} \cdot 1.034
\end{align*}
where $H,L,C$ and $O$ indicate, respectively, the high, low, close and open BTC-USD rate for the day. In our application, the dependent variable is the natural logarithm of the Garman-Klass proxy.  Figure \ref{fig:series} shows its evolution during the analysis period, the Normal QQ plot and the (partial) autocorrelation function.

\begin{figure}[h!]
\centering
\caption{Garman-Klass volatility proxy in log scale. Panel (A) displays the time series evolution during the analysis period; Panel (B) shows the evolution of the time series starting from June 2017 and highlights the announcement dates and the launch of the two futures (see Table \ref{tab:event_history} for the details and the exact dates); Panel (c) shows the Normal QQ Plot, autocorrelation and partial autocorrelation functions.}
\label{fig:series}
\begin{tabular}{m{1cm}m{15cm}}
A) & \includegraphics[scale=0.35]{series_all} \\
B) & \includegraphics[scale=0.35]{series} \\
C) & \includegraphics[scale=0.35]{logGK_plots} \\
\end{tabular}
\end{figure}

In order to shed light on the driving forces underneath Bitcoin volatility, we also investigate the impact of Bitcoin futures on the daily transaction volumes. Indeed, the possible volatility surge (decline) following the launch of the two contracts might be driven by increased (decreased) volumes. From Figure \ref{fig:volumes}, showing the evolution of daily volumes throughout the analysis period, we can notice that Bitcoin transactions experienced a sharp increase during 2014. However, as shown in Panel (B), by focusing on a restricted time period, the time series becomes more stable.

\begin{figure}[h!]
\centering
\caption{Bitcoin daily transaction volumes in log scale. Panel (A) displays the time series evolution during the analysis period; Panel (B) shows the evolution of the time series starting from February 2015 and highlights the announcement dates and the launch of the two futures (see Table \ref{tab:event_history} for the details and the exact dates); Panel (c) shows the Normal QQ Plot, autocorrelation and partial autocorrelation functions.}
\label{fig:volumes}
\begin{tabular}{m{1cm}m{15cm}}
A) & \includegraphics[scale=0.35]{volumes_1} \\
B) & \includegraphics[scale=0.35]{volumes_2} \\
C) & \includegraphics[scale=0.35]{volumes_plot} \\
\end{tabular}
\end{figure}

\clearpage
\subsection{Methodology}
\label{subsect:method}
As described in Section \ref{subsect:estimators}, to estimate the effect on Bitcoin volatility and volumes of the two futures and each related announcement, we need to fit a C-ARIMA model up to the day before each intervention, which, in our case, leads to the estimation of  six different models: based on the first four, we make $1$-step ahead predictions to compute the contemporaneous effect of each announcement, whereas, based on the last two models, we compute the effect of futures at different time horizons. In particular, since the CBOE future was launched $1$ week before the CME future, for the former we estimate the temporal average pointwise effect at $1$-week horizon, whereas for the CME futures we can also consider the $2$-week and the $3$-week horizons.

Notice that since both dependent variables are in log scale, the difference between the observed and the predicted outcomes is equal to the log of their ratio. This means that we are assuming a multiplicative effect and, as a result, among the estimators defined by Equations (\ref{eqn:point_estimator})--(\ref{eqn:avgpoint_estimator}) we can focus on the temporal average pointwise effect. Indeed, it has a financial interpretation after re-exponentiating, being the geometric average of the point effects.

Volatility is typically stationary and non-seasonal, hence, our six independent models are all built from an ARMA$(p,q)$. The models also include covariates, so to improve the forecast of Bitcoin volatility in the absence of intervention. We selected the set of covariates based on a survey of relevant literature (see Section \ref{sect:background}). In particular, we included the daily log returns and the Garman-Klass proxies of the EUR-USD exchange rate, the MSCI Emerging Market Index, the Shanghai Stock Exchange Composite Index and the MSCI World Index. Then, among the economic factors we selected the Federal Reserve money supply M1 aggregate, the US monthly inflation rate, the Federal funds Target Rate and the US GDP. Technology factors are represented by the hash rate (estimated number of tera hashes per second the Bitcoin network is performing) and Bitcoin supply (total number of Bitcoins in circulation). All covariates are in log scale and have been made stationary in case they are not\footnote{First differences have been taken for \texttt{m1}, \texttt{midrate}, \texttt{gdp} and \texttt{hash}, whereas \texttt{btcsupply} has been differenced twice.}. Finally, since the second Bitcoin halving took place on July 9, 2016, we also included a dummy variable taking value $1$ after the halving. Table \ref{tab:corr} summarizes the correlation between the covariates and the Garman-Klass volatility. Overall, the Garman-Klass proxy seems to have a small correlation with the volatility of the EUR-USD exchange rate. All the considered covariates are reasonably unaffected by Bitcoin futures; thus, we can use their values post-intervention to improve the prediction accuracy of the counterfactual outcomes absent the futures.

The selection of the six independent models was based on the Bayesian Information Criterion (BIC). Since the announcement dates are close to each other and to the actual launch of the futures, the characteristics of the data are in all cases well described by the same C-ARMA$(2,1)$ specification, i.e.,\footnote{Notice from Figure \ref{fig:series} that the acf of the Bitcoin Garman-Klass volatility indicator shows a very slow decay, a behavior which is sometimes labeled as \textit{long-memory}.
To capture this pattern, specific models have been developed, like the HAR (Heterogeneus Autoregressive) and its log-counterpart, the log-HAR \citep{Corsi:2009}. On the other hand, \citet[][Sect.~5.3]{Cipollini:Gallo:Palandri:2020} demonstrate that some (2,1)-specifications of different ARMA-like models replicate the ability to approximate the long memory pattern observed in the autocorrelation of realized variances, a feature which has made the HAR model popular.
Such an ability is due to the high estimated \textit{persistence} and the presence of a second order parameter, usually significantly negative, in the AR part.}
\begin{equation}
\label{eqn:volatility}
\Y_t(\w) = \frac{\theta_1(L)}{\phi_2(L)}\varepsilon_t + \X_t' \beta + \sum_{j = 1}^{n-1} \tau_t^{(j)} + \tau_t^{(n)}\mathds{1}_{\{\w=1 \}} .
\end{equation}
Unlike the Garman-Klass proxy, log-volumes are not stationary (the KPSS test rejects the null hypothesis of stationarity at the $1\%$ level) and show a weekly seasonal pattern; thus, denoting with $\V_t(\w)$ the potential outcome time series of log-volumes at time $t$, the model selected with the BIC criterion is a C-ARIMA$(1,1,1)(0,0,2)_7$, which is written as,
\begin{equation}
\label{eqn:volumes}
\V_t(\w) = \frac{\Theta_2(L^7) \theta_1(L)}{(1-L)\phi_1(L)}\varepsilon_{t} + \X_t' \beta + \sum_{j = 1}^{n-1} \tau_t^{(j)} + \tau_t^{(n)}\mathds{1}_{\{\w=1 \}}.
\end{equation}
In case that trading volumes turn out to be affected by the two futures and if the volume-volatility relationship holds also in the Bitcoin market, the effect on Bitcoin volatility estimated by (\ref{eqn:volatility}) would then be due, at least in part, by increased (decreased) transactions.
To find the impact on Bitcoin volatility attributable solely to the two futures, we should be able to predict the counterfactual outcome that we would have had absent the futures and their impact on volumes. We can do that by using the following new model for the Garman-Klass volatility proxy:
\begin{equation}
\label{eqn:mixed_model}
\Y_t(\w) = \frac{\theta_q(L)}{\phi_p(L)}\varepsilon_{t} + \X_t' \beta + (1-L) \sum_{i = 1}^3 \alpha_i \tilde{\V}_{t-i} + \sum_{j = 1}^{n-1} \tau_t^{(j)} + \tau_t^{(n)}\mathds{1}_{\{\w=1 \}}
\end{equation}
where $\tilde{\V}_{t} = (\V_1, \dots, \V_{t_n-1},\hat{\V}_{t_n},\dots,\hat{\V}_T)$ is a vector formed by the observed volumes up to the day before the $n$-th intervention and then by the counterfactual volumes forecasted by model (\ref{eqn:volumes}), denoted with $\hat{\V}_t$; thus, $\tilde{\V}_{t-1}$, $\tilde{\V}_{t-2}$ and $\tilde{\V}_{t-3}$ are the lagged values of $\tilde{\V}_{t}$.\footnote{Notice that including $\tilde{\V}_t$ in the set of covariates does not violate the covariates-treatment independence assumption. Indeed, by using the predicted counterfactual volumes $\hat{\V}_t$ in the post intervention period, $\tilde{\V}_t$ is truly unaffected by the interventions. } The idea is that, in the pre-intervention period, significant coefficients in the lagged volumes would confirm the presence of a volume-volatility relationship in the Bitcoin market, whereas, in the post-intervention period, the lagged predicted volumes would improve the estimation of the counterfactual volatility absent the futures and their effect on volumes. The next Section reports the results of the empirical analysis on the Garman-Klass volatility proxy and on Bitcoin daily volumes.

\section{Results}
\label{sect:results}
The parameter estimates of the six C-ARMA$(2,1)$ models, the estimated contemporaneous effect of the announcements and the estimated temporal average pointwise effects of the two futures on the Garman-Klass volatility proxy are reported in Table \ref{tab:estimates}. Since the residuals diagnostics do not seem to support the Normality assumption (see Figure \ref{fig:residuals_GK}), we report empirical critical values computed from Equation (\ref{eqn:avgpoint_estimator_var_h0}) by bootstrapping the error terms from the model residuals.\footnote{The number of bootstrap samples is $10,000$. See Table \ref{tab:sd_norm} for the estimated standard errors of the temporal average pointwise effects computed from Equation (\ref{eqn:avgpoint_estimator_var}).}
We find no evidence of significant contemporaneous effects of the announcements, meaning that the news of upcoming futures was not sufficient to spark investors' interests in the (crypto)currency. Instead, we observe an interesting result for the actual introduction of the two futures: there is no evidence of a significant effect related to CBOE futures, whereas the effects related to CME futures are significant at all time horizons. For example, we find that, $1$-week after its launch, the CME future increased volatility by, on average, $+116\%$.
The analysis on Bitcoin volumes might provide an explanation of the different results found for the CBOE and CME instruments.
Figure \ref{fig:pointwise_cme} shows the forecasted series and the pointwise causal effects. Having found no evidence of association between the covariates and the Garman-Klass proxy (lagged values of the volatility seem sufficient in describing the dynamics of current Bitcoin volatility), as a robustness check to the discussed results we also estimated six alternative models with no regressors. The estimated causal effects based on the alternative models are in line with the results reported in this section (see Table \ref{tab:estimates_2} for the details).

Table \ref{tab:estimates_volumes_log} shows the results of the analysis performed on Bitcoin daily transaction volumes (in log scale). Again, we report empirical critical values computed by bootstrapping the errors.\footnote{The number of bootstrap samples is $10,000$. See Table \ref{tab:sd_norm} for the estimated standard errors of the temporal average pointwise effects computed from Equation (\ref{eqn:avgpoint_estimator_var}).} This time we found significant effects for both the CBOE and the CME futures at all time horizons, whereas there is no evidence of a contemporaneous causal effect of the announcements. Trading volumes show a positive association with Bitcoin supply, meaning that the growing number of coins also translates into more transactions. In addition, Bitcoin volumes seem to be negatively related to the MSCI World Index, suggesting that the cryptocurrency might act as a safe heaven during turbulence in traditional markets. Figure \ref{fig:pointwise_cme_vol} provides a graphical representation of these results. As a robustness check, we also repeated the analyses on a shorter time interval starting from February 2015 when Bitcoin volumes, after a sharp increase during 2014, finally reached a more stable path. The results are in line with those reported in this section and are shown in Table \ref{tab:estimates_volumes_2}. Interestingly, the effects of the two futures have opposite sign: transaction volumes decreased of approximately $-18\%$ due to the launch of the CBOE contract, which might indicate that for a short period of time investors privileged the derivative over of the underlying. Lower transaction volumes might also explain the absence of effect on volatility. Conversely, the CME contract increased volumes by $+39\%$ in the first week after its launch, suggesting that, benefiting from the increased transparency of the market, investors' interest toward Bitcoin rose (especially in the short-run) boosting Bitcoin trading and, in turn, its volatility.

To check whether the above interpretation is correct and to find the impact on volatility due solely to the futures, we estimated one last model for the Garman-Klass volatility proxy where we included a few lagged values of transaction volumes in the set of predictors, as in Equation (\ref{eqn:mixed_model}): if the effect on volatility was, at least in part, triggered by volumes, we would expect to find significant coefficients for the lagged log-volumes (since volumes exhibit non-stationarity, we considered their first difference). Then, as described in Section \ref{subsect:method}, to improve the prediction of the counterfactual volatility absent the futures and their impact on volumes, in the post-intervention period we used $\hat{\V_t}$. Table \ref{tab:mediation} reports the results: as expected, the coefficient of lagged volumes is positive and significant, suggesting that when transactions increase Bitcoin volatility rises as well. Interestingly, past volumes seem to exert a strong influence on volatility even at $t-3$.\footnote{We also tested the presence of possible feedback effects stemming from volatility to volumes. To do that, in Equation (\ref{eqn:mixed_model}) we switched lagged log-volumes with lagged log-volatilities and $\Y_t$ with $\V_t$. As shown in Table \ref{tab:mediation_volatility}, we did not find evidence of a feedback effect.} Notice that by considering the volumes that we would have observed in the absence of intervention, the estimated causal effect on Bitcoin volatility reported in Table \ref{tab:mediation} is the actual impact due solely to the two futures, which is then higher than the effect reported in Table \ref{tab:estimates}. For example, we obtain that, considering the lower volumes that we would have had in the absence of the future, the CME contract is actually responsible of a $+139\%$ increase in volatility in first week of its launch
(instead of the $+116\%$ estimated by the first model).\footnote{This is because, by considering lower transactions, the counterfactual volatility predicted by model (\ref{eqn:mixed_model}) is below the one predicted by model (\ref{eqn:volatility}).}  Again, the effect of the CBOE contract was not significant, whereas the causal effect of the CME contract is significant at all time horizons, with a small reduction after $3$ weeks.

\begin{table}[h!]
\centering
\caption{Estimates of the six independent C-ARMA models fitted on the historical daily values of the \textbf{Garman-Klass volatility proxy} (in log scale) starting from May 3, 2014 and up to the day before each intervention (standard errors within parentheses). See Table \ref{tab:event_history} for the exact dates of the announcements and futures launches. In this table, $\hat{\tau}^{(n)}$ indicates the estimated contemporaneous effect of each announcement; $\hat{\tau}^{(CBOE)}$ is the temporal average pointwise effect of the CBOE futures; and, $\hat{\tau}_t^{(CME)}$ is the temporal average pointwise effect of the CME futures at $1$-week, $2$-weeks and $3$-weeks horizons, respectively $t = 7$, $t = 14$ and $t = 21$. The multiplicative effects can be recovered by re-exponentiating the estimated effects. The empirical critical values for the causal effects are computed from Equation (\ref{eqn:avgpoint_estimator_var_h0}) by bootstrapping the errors from the residuals of each model.}
\label{tab:estimates}
\scalebox{0.75}{
\input{Tables/BTC_gk_est}
}
\end{table}

\begin{table}[ht]
\centering
\caption{Estimates of six C-ARIMA models fitted on the historical \textbf{Bitcoin daily volumes} (in log scale) starting from May 3, 2014 and up to the day before each intervention (standard errors within parentheses). See Table \ref{tab:event_history} for the exact dates of the announcements and futures launches. In this table, $\hat{\tau}^{(n)}$ indicates the estimated contemporaneous effect of each announcement; $\hat{\tau}^{(CBOE)}$ is the temporal average pointwise effect of the CBOE futures; and, $\hat{\tau}_t^{(CME)}$ is the temporal average pointwise effect of the CME futures at $1$-week, $2$-weeks and $3$-weeks horizons, respectively $t = 7$, $t = 14$ and $t = 21$). The multiplicative effects can be recovered by re-exponentiating the estimated effects. The empirical critical values for the causal effects are computed from Equation (\ref{eqn:avgpoint_estimator_var_h0}) by bootstrapping the errors from the residuals of each model.}
\label{tab:estimates_volumes_log}
\scalebox{0.75}{
\input{Tables/BTC_volumes_est}
}
\end{table}

\begin{table}
\centering
\caption{Estimates of the C-ARMA models fitted on the historical daily values of the \textbf{Garman-Klass volatility proxy} (in log scale) with lagged volumes among regressors (standard errors within parentheses). The models are estimated starting from May 3, 2014 and up to the day before each intervention. See Table \ref{tab:event_history} for the exact dates of the futures launches. In this table, $\hat{\tau}^{(CBOE)}$ is the temporal average pointwise effect of the CBOE futures and $\hat{\tau}_t^{(CME)}$ is the temporal average pointwise effect of the CME futures at $1$-week, $2$-weeks and $3$-weeks horizons (indicated with $t = 7$, $t = 14$ and $t = 21$). The multiplicative effects can be recovered by re-exponentiating the estimated effects. The empirical critical values for the causal effects are computed from Equation (\ref{eqn:avgpoint_estimator_var_h0}) by bootstrapping the errors from the residuals of each model.}
\label{tab:mediation}
\scalebox{0.75}{
\input{Tables/BTC_mediation}
}
\end{table}

\newgeometry{bottom=1.5cm,top=1.5cm,left=1cm,right=1cm}
\begin{figure}[h!]
\centering
\caption{Observed and forecasted \textbf{Garman-Klass volatility proxy} (in log scale) at $1$-week, $2$-weeks and $3$-weeks from the launch of the CME future, indicated by the dashed bar. The right charts show the resulting effects (computed as the difference between the observed and the forecasted series) with their $95 \%$ confidence bounds.}
\label{fig:pointwise_cme}
\begin{tabular}{cc}
\includegraphics[scale=0.40]{forecast_GK} & \includegraphics[scale=0.40]{impact_GK} \\
\end{tabular}
\end{figure}

\begin{figure}[h!]
\centering
\caption{Observed and forecasted \textbf{Bitcoin daily volumes} (in log scale) at $1$-week, $2$-weeks and $3$-weeks from the launch of the CME future, indicated by the dashed bar. The right charts show the resulting effects (computed as the difference between the observed and the forecasted series) with their $95 \%$ confidence bounds.}
\label{fig:pointwise_cme_vol}
\begin{tabular}{cc}
\includegraphics[scale=0.4]{forecast_volumes} & \includegraphics[scale=0.4]{impact_volumes} \\
\end{tabular}
\end{figure}
\restoregeometry

\newgeometry{bottom=2.5cm,top=2.3cm,left=2cm,right=2cm}
\section{Conclusions}
\label{sect:conclusions}
In December 2017, two leading derivative exchanges, CBOE and CME, introduced the first two regulated futures having Bitcoin as underlying asset. Contributing to the nascent stream of literature seeking to assess the effect generated by the new instruments, the goal of this paper was estimating their causal impact on Bitcoin volatility; we also investigated the presence of announcement effects as well as the impact on Bitcoin trading volumes.

To estimate the causal effect generated by the launch of Bitcoin futures,
we employed a novel methodology, C-ARIMA, based on the Rubin's potential outcomes framework. After a detailed discussion of the assumptions enabling the definition of the effect and its attribution to the a specific treatment, this approach infers the causal effect of an intervention by direct comparison between the observed outcome and a predicted counterfactual. In this paper, we extended the C-ARIMA approach to a multi-intervention setting by formalizing an additional assumption and introducing two causal estimands capturing the effect of the multiple announcements and of the futures introductions. Then, we defined estimators of these effects under the C-ARIMA model and derived hypothesis tests.

The results indicate that the CME contract produced an increase of both Bitcoin volatility and trading volume. We also found evidence that lagged volumes are positively associated to Bitcoin volatility, suggesting that the effect of the introduction of CME futures on volatility might be explained, at least partially, by the increased volumes following the same intervention. After controlling for the lagged predicted volumes in the absence of intervention, we were able to find the effect on volatility actually attributable to the future contract. This effect was positive and significant at all time horizons. Conversely, the CBOE futures had no impact on volatility and produced a small negative effect on volumes. Finally, we did not find evidence of announcement effects.

Overall, our results support the general finding in the literature of a positive association between volumes and volatility. Furthermore, they show that future trading (although in a regulated environment) contributed to increase Bitcoin volatility instead of diminishing it. Whether this is the result of improved transparency in the market (promoting transactions and hence producing more volatile returns) or increased speculation, this important evidence is backed up by a reliable causal approach: regulators should then be aware of these mechanisms as banks and financial institutions get increasingly involved in the crypto world.

\restoregeometry

\newpage