{
"cells": [
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"# Synthetic control with `pymc` models"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import arviz as az\n",
"\n",
"import causalpy as cp"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"%load_ext autoreload\n",
"%autoreload 2\n",
"%config InlineBackend.figure_format = 'retina'\n",
"seed = 42"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"## Load data"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"df = cp.load_data(\"sc\")\n",
"treatment_time = 70"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"## Run the analysis\n",
"\n",
":::{note}\n",
"The `random_seed` keyword argument for the PyMC sampler is not necessary. We use it here so that the results are reproducible.\n",
":::"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" a \n",
" b \n",
" c \n",
" d \n",
" e \n",
" f \n",
" g \n",
" counterfactual \n",
" causal effect \n",
" actual \n",
" \n",
" \n",
" \n",
" \n",
" 0 \n",
" 0.793234 \n",
" 1.277264 \n",
" -0.055407 \n",
" -0.791535 \n",
" 1.075170 \n",
" 0.817384 \n",
" -2.607528 \n",
" 0.144888 \n",
" -0.0 \n",
" 0.398287 \n",
" \n",
" \n",
" 1 \n",
" 1.841898 \n",
" 1.185068 \n",
" -0.221424 \n",
" -1.430772 \n",
" 1.078303 \n",
" 0.890110 \n",
" -3.108099 \n",
" 0.601862 \n",
" -0.0 \n",
" 0.491644 \n",
" \n",
" \n",
" 2 \n",
" 2.867102 \n",
" 1.922957 \n",
" -0.153303 \n",
" -1.429027 \n",
" 1.432057 \n",
" 1.455499 \n",
" -3.149104 \n",
" 1.060285 \n",
" -0.0 \n",
" 1.232330 \n",
" \n",
" \n",
" 3 \n",
" 2.816255 \n",
" 2.424558 \n",
" 0.252894 \n",
" -1.260527 \n",
" 1.938960 \n",
" 2.088586 \n",
" -3.563201 \n",
" 1.520801 \n",
" -0.0 \n",
" 1.672995 \n",
" \n",
" \n",
" 4 \n",
" 3.865208 \n",
" 2.358650 \n",
" 0.311572 \n",
" -2.393438 \n",
" 1.977716 \n",
" 2.752152 \n",
" -3.515991 \n",
" 1.983661 \n",
" -0.0 \n",
" 1.775940 \n",
" \n",
" \n",
"
\n",
"
"
]
},
"metadata": {
"image/png": {
"height": 811,
"width": 711
}
},
"output_type": "display_data"
}
],
"source": [
"fig, ax = result.plot(plot_predictors=True)"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"================================SyntheticControl================================\n",
"Control units: ['a', 'b', 'c', 'd', 'e', 'f', 'g']\n",
"Treated unit: actual\n",
"Model coefficients:\n",
" a 0.34, 94% HDI [0.29, 0.38]\n",
" b 0.049, 94% HDI [0.0087, 0.09]\n",
" c 0.31, 94% HDI [0.26, 0.35]\n",
" d 0.055, 94% HDI [0.01, 0.1]\n",
" e 0.025, 94% HDI [0.0011, 0.067]\n",
" f 0.19, 94% HDI [0.11, 0.26]\n",
" g 0.038, 94% HDI [0.0025, 0.088]\n",
" sigma 0.26, 94% HDI [0.22, 0.31]\n"
]
}
],
"source": [
"result.summary()"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"As well as the model coefficients, we might be interested in the average causal impact and average cumulative causal impact."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" mean \n",
" sd \n",
" hdi_3% \n",
" hdi_97% \n",
" mcse_mean \n",
" mcse_sd \n",
" ess_bulk \n",
" ess_tail \n",
" r_hat \n",
" \n",
" \n",
" \n",
" \n",
" x[actual] \n",
" -1.714 \n",
" 0.217 \n",
" -2.141 \n",
" -1.332 \n",
" 0.007 \n",
" 0.004 \n",
" 901.0 \n",
" 1428.0 \n",
" 1.0 \n",
" \n",
" \n",
"
\n",
"
\n",
"\n",
"
\n",
" \n",
" \n",
" mean \n",
" sd \n",
" hdi_3% \n",
" hdi_97% \n",
" mcse_mean \n",
" mcse_sd \n",
" ess_bulk \n",
" ess_tail \n",
" r_hat \n",
" \n",
" \n",
" \n",
" \n",
" x[actual] \n",
" -51.424 \n",
" 6.525 \n",
" -64.234 \n",
" -39.965 \n",
" 0.215 \n",
" 0.109 \n",
" 901.0 \n",
" 1428.0 \n",
" 1.0 \n",
" \n",
" \n",
"
\n",
"