{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "3d83155f",
   "metadata": {
    "kernel": "Stata"
   },
   "source": [
    "# Simulating and Estimating Spillovers in Stata\n",
    "\n",
    "This notebook simulates and then estimates spillover effects where potential outcomes are described by the function\n",
    "$$\n",
    "y(\\tau_{i,t}, v(\\tau_{-i,t}))\n",
    "$$\n",
    "\n",
    "for firm-level treatment $\\tau_{i,t}$.  Outcomes depend on peer firm treatment status through the scalar function $v(\\tau_{-i,t})$ where $\\tau_{-i,t}$ denotes treatments outside of the focal firm.  The scalar function was first devised by Hong and Raudenbush (2006 JASA) to address the dimensionality of $\\tau_{-i,1}$ (the function $y(\\tau_{i,t},\\tau_{-i,t})$ has $2^N-1$ potential outcomes).  The scalar function $v(\\cdot)$ is specified as\n",
    "$$\n",
    "v(\\tau_{-i,t}) := \\frac{1}{|C|}\\sum_{i' \\in C} \\tau_{i',t}\n",
    "$$\n",
    "for $C := \\{i'|\\text{Ind. }i'=\\text{Ind. }i\\}$.  This follows Ferraci, Jolivet, and van den Berg (2014 REStat).  As shorthand, we will write $v(\\tau_{-i,t})$ as $\\rho_{j,t}$ where $j$ indexes the industry.\n",
    "\n",
    "We begin by constructing a panel with $40$ industries, $50$ firms per industry, and $40$ quarters of data.  A handful of true $\\beta$ parameters are also specified."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "150217d8",
   "metadata": {
    "kernel": "Stata"
   },
   "outputs": [],
   "source": [
    "// clear everything from Stata's memory\n",
    "clear\n",
    "\n",
    "// create industries\n",
    "qui set obs 40\n",
    "gen j = _n\n",
    "\n",
    "// for each industry, create a set of firms\n",
    "qui expand 50\n",
    "bys j: gen i = j*100+_n\n",
    "\n",
    "// for each industry-firm, create a set of time periods\n",
    "qui expand 40\n",
    "bys j i: gen t = _n\n",
    "\n",
    "// true model parameters\n",
    "scalar b0 = 0\n",
    "scalar b1 = -3\n",
    "scalar b2 = 2\n",
    "scalar b3 = -2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f2a814c5",
   "metadata": {
    "kernel": "Stata"
   },
   "source": [
    "### Exogenous Treatment Example\n",
    "\n",
    "Suppose first that treatment assignment is random and that spillover effects are linear. This implies a model\n",
    "$$\n",
    "y_{i,t} = \\beta_0 + \\beta_1\\tau_{i,t} + \\beta_2\\tau_{i,t}\\rho_{j,t} + \\beta_3(1-\\tau_{i,t})\\rho_{j,t} + \\epsilon_{i,t}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "9aba668a",
   "metadata": {
    "kernel": "Stata"
   },
   "outputs": [],
   "source": [
    "// reproducibility\n",
    "set seed 0\n",
    "\n",
    "// create firm and industry treatment\n",
    "gen tau = rbinomial(1, .05)\n",
    "bys j t: egen rho = mean(tau)\n",
    "\n",
    "// create a true model\n",
    "gen u = rnormal()\n",
    "gen y = b0 + b1*tau + b2*1.tau*rho + b3*0.tau*rho + u"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3bdc1929",
   "metadata": {
    "kernel": "Stata"
   },
   "source": [
    "When treatment is exogenous and spillover effects are linear in $\\rho$, treatment effects are easy to estimate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "2c36bd25",
   "metadata": {
    "kernel": "Stata"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "      Source |       SS           df       MS      Number of obs   =    80,000\n",
      "-------------+----------------------------------   F(3, 79996)     =   9892.12\n",
      "       Model |  29701.2931         3  9900.43102   Prob > F        =    0.0000\n",
      "    Residual |  80063.2094    79,996  1.00084016   R-squared       =    0.2706\n",
      "-------------+----------------------------------   Adj R-squared   =    0.2706\n",
      "       Total |  109764.502    79,999  1.37207343   Root MSE        =    1.0004\n",
      "\n",
      "------------------------------------------------------------------------------\n",
      "           y | Coefficient  Std. err.      t    P>|t|     [95% conf. interval]\n",
      "-------------+----------------------------------------------------------------\n",
      "         tau |  -2.952744   .0389518   -75.80   0.000    -3.029089   -2.876398\n",
      "             |\n",
      "   tau#c.rho |\n",
      "          0  |  -2.218683   .1155108   -19.21   0.000    -2.445084   -1.992283\n",
      "          1  |   .8730354   .4987173     1.75   0.080    -.1044473    1.850518\n",
      "             |\n",
      "       _cons |   .0107455   .0066978     1.60   0.109    -.0023821     .023873\n",
      "------------------------------------------------------------------------------\n"
     ]
    }
   ],
   "source": [
    "reg y tau i.tau#c.rho"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "56907b6d",
   "metadata": {
    "kernel": "Stata"
   },
   "source": [
    "The Stable Unit Treatment Value Assumption (SUTVA) imposes the condition that $y(\\tau_{i,t},\\tau_{-i,t})=y(\\tau_{i,t})$ so that there are no spillover effects.  This implies the model\n",
    "$$\n",
    "y_{i,t} = \\theta_0 + \\theta_1\\tau_{i,t} + u_{i,t}\n",
    "$$\n",
    "\n",
    "Theoretically, $\\hat{\\theta_1} = \\hat{\\beta_1} + (\\hat{\\beta_2}-\\hat{\\beta_3})\\bar{\\rho}_{j,t}$.  Within the simulated sample, the observed bias is close to the asymptotic bias."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "e3ad1066",
   "metadata": {
    "kernel": "Stata"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "\n",
      "      Source |       SS           df       MS      Number of obs   =    80,000\n",
      "-------------+----------------------------------   F(1, 79998)     =  29169.45\n",
      "       Model |  29328.9848         1  29328.9848   Prob > F        =    0.0000\n",
      "    Residual |  80435.5176    79,998  1.00546911   R-squared       =    0.2672\n",
      "-------------+----------------------------------   Adj R-squared   =    0.2672\n",
      "       Total |  109764.502    79,999  1.37207343   Root MSE        =    1.0027\n",
      "\n",
      "------------------------------------------------------------------------------\n",
      "           y | Coefficient  Std. err.      t    P>|t|     [95% conf. interval]\n",
      "-------------+----------------------------------------------------------------\n",
      "         tau |  -2.783437   .0162974  -170.79   0.000    -2.815379   -2.751494\n",
      "       _cons |  -.0973879   .0036369   -26.78   0.000    -.1045162   -.0902596\n",
      "------------------------------------------------------------------------------\n",
      "\n",
      "\n",
      "Theoretical bias .1991999967209996\n",
      "\n",
      "Observed bias .2165633739491448\n"
     ]
    }
   ],
   "source": [
    "reg y tau\n",
    "\n",
    "qui su rho\n",
    "di \"Theoretical bias `=(b2-b3)*`r(mean)''\"\n",
    "di \"Observed bias `=_b[tau]-b1'\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8e5c6eb7",
   "metadata": {
    "kernel": "Stata"
   },
   "source": [
    "### Endogenous Treatment Sample\n",
    "\n",
    "Now suppose that treatment assignment is endogenous.  In particular, there is an industry-level covariate $w$ with which a firm-level covariate $x$ is correlated.  We allow the firm-level covariate $x$ to predict treatment status.\n",
    "\n",
    "For simplicity, we won't allow $x$ to affect outcomes directly, though adding this feature wouldn't really change anything."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "dd857775",
   "metadata": {
    "kernel": "Stata"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "\n",
      "\n",
      "\n",
      "\n",
      "\n",
      "\n",
      "(39,503 missing values generated)\n",
      "\n",
      "(39,503 missing values generated)\n",
      "\n",
      "\n",
      "\n",
      "\n",
      "\n",
      "(39,503 missing values generated)\n",
      "\n",
      "(39,503 missing values generated)\n"
     ]
    }
   ],
   "source": [
    "// reproducibility\n",
    "set seed 0\n",
    "\n",
    "// hold on only to the panel identifiers\n",
    "keep i j t\n",
    "\n",
    "// create dummy for unique industry-quarters\n",
    "bys j t: gen _first = _n == 1\n",
    "\n",
    "// generate industry-quarter data\n",
    "qui gen w = runiform(1, 3) if _first == 1\n",
    "sort j t i\n",
    "qui by j t: carryforward w, replace\n",
    "\n",
    "// generate firm data\n",
    "gen x = rnormal(w, 1)\n",
    "\n",
    "// simulate treatment\n",
    "gen pr = 1/(1+exp(2*x))\n",
    "\n",
    "// create treatment dummy\n",
    "qui gen tau = rbinomial(1,pr)\n",
    "\n",
    "// rbinomial yields some missing values for pr very close to zero or one\n",
    "qui replace tau = 0 if mi(tau) & pr < .5\n",
    "qui replace tau = 1 if mi(tau) & pr >= .5\n",
    "\n",
    "// industry treatment intensity\n",
    "bys j t: egen rho = mean(tau)\n",
    "\n",
    "// generate outcome of interest\n",
    "gen u = rnormal() + .05*x\n",
    "gen y = b0 + b1*tau + b2*1.tau*rho + b3*0.tau*rho + u"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "47419699",
   "metadata": {
    "kernel": "Stata"
   },
   "source": [
    "As one would expect, without addressing endogenous treatment assignment, the model that we ran before is flawed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "a9921337",
   "metadata": {
    "kernel": "Stata"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "      Source |       SS           df       MS      Number of obs   =    40,497\n",
      "-------------+----------------------------------   F(3, 40493)     =   4911.02\n",
      "       Model |  14675.0009         3  4891.66695   Prob > F        =    0.0000\n",
      "    Residual |  40333.4624    40,493  .996060118   R-squared       =    0.2668\n",
      "-------------+----------------------------------   Adj R-squared   =    0.2667\n",
      "       Total |  55008.4632    40,496  1.35836782   Root MSE        =    .99803\n",
      "\n",
      "------------------------------------------------------------------------------\n",
      "           y | Coefficient  Std. err.      t    P>|t|     [95% conf. interval]\n",
      "-------------+----------------------------------------------------------------\n",
      "         tau |  -3.071996   .0404585   -75.93   0.000    -3.151296   -2.992697\n",
      "             |\n",
      "   tau#c.rho |\n",
      "          0  |   -2.02764   .0225631   -89.87   0.000    -2.071864   -1.983415\n",
      "          1  |   2.036824   .0798152    25.52   0.000     1.880384    2.193263\n",
      "             |\n",
      "       _cons |   .1169907   .0099473    11.76   0.000     .0974938    .1364876\n",
      "------------------------------------------------------------------------------\n"
     ]
    }
   ],
   "source": [
    "// without addressing endogeneity\n",
    "reg y tau i.tau#c.rho"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ba734360",
   "metadata": {
    "kernel": "Stata"
   },
   "source": [
    "The simplest approach is to use regression adjustment and include $x$ as a covariate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "31987b84",
   "metadata": {
    "kernel": "Stata"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "      Source |       SS           df       MS      Number of obs   =    40,497\n",
      "-------------+----------------------------------   F(4, 40492)     =   3713.66\n",
      "       Model |  14763.8806         4  3690.97015   Prob > F        =    0.0000\n",
      "    Residual |  40244.5826    40,492  .993889722   R-squared       =    0.2684\n",
      "-------------+----------------------------------   Adj R-squared   =    0.2683\n",
      "       Total |  55008.4632    40,496  1.35836782   Root MSE        =    .99694\n",
      "\n",
      "------------------------------------------------------------------------------\n",
      "           y | Coefficient  Std. err.      t    P>|t|     [95% conf. interval]\n",
      "-------------+----------------------------------------------------------------\n",
      "         tau |   -2.99079   .0413167   -72.39   0.000    -3.071772   -2.909809\n",
      "             |\n",
      "   tau#c.rho |\n",
      "          0  |  -2.006762   .0226464   -88.61   0.000    -2.051149   -1.962374\n",
      "          1  |   2.053711   .0797482    25.75   0.000     1.897403     2.21002\n",
      "             |\n",
      "           x |   .0452351   .0047835     9.46   0.000     .0358594    .0546108\n",
      "       _cons |   .0105126   .0150171     0.70   0.484    -.0189213    .0399465\n",
      "------------------------------------------------------------------------------\n"
     ]
    }
   ],
   "source": [
    "// regression adjustment\n",
    "reg y tau i.tau#c.rho x"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2babb097",
   "metadata": {
    "kernel": "Stata"
   },
   "source": [
    "Now suppose that spillover effects are nonlinear in $\\rho$.\n",
    "\n",
    "In this scenario, regression adjustment alone isn't sufficient."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "6139e965",
   "metadata": {
    "kernel": "Stata"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "(40,497 real changes made)\n",
      "\n",
      "\n",
      "      Source |       SS           df       MS      Number of obs   =    40,497\n",
      "-------------+----------------------------------   F(4, 40492)     =  36948.53\n",
      "       Model |  159000.108         4  39750.0271   Prob > F        =    0.0000\n",
      "    Residual |  43562.1663    40,492  1.07582155   R-squared       =    0.7849\n",
      "-------------+----------------------------------   Adj R-squared   =    0.7849\n",
      "       Total |  202562.275    40,496  5.00203167   Root MSE        =    1.0372\n",
      "\n",
      "------------------------------------------------------------------------------\n",
      "           y | Coefficient  Std. err.      t    P>|t|     [95% conf. interval]\n",
      "-------------+----------------------------------------------------------------\n",
      "         tau |  -9.582966   .0429859  -222.93   0.000    -9.667219   -9.498712\n",
      "             |\n",
      "   tau#c.rho |\n",
      "          0  |  -4.528388   .0235614  -192.20   0.000    -4.574568   -4.482207\n",
      "          1  |   4.051024   .0829701    48.83   0.000     3.888401    4.213648\n",
      "             |\n",
      "           x |   .0640858   .0049767    12.88   0.000     .0543313    .0738403\n",
      "       _cons |   3.418386   .0156238   218.79   0.000     3.387763     3.44901\n",
      "------------------------------------------------------------------------------\n"
     ]
    }
   ],
   "source": [
    "// nonlinear spillovers\n",
    "replace y = b0 + b1*tau + b2*1.tau*log(.1+rho) + b3*0.tau*log(.1+rho) + u\n",
    "\n",
    "// regression adjustment\n",
    "reg y tau i.tau#c.rho x"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1bf45019",
   "metadata": {
    "kernel": "Stata"
   },
   "source": [
    "Begin by computing likelihoods of treatment.\n",
    "\n",
    "Ideally, this would be done on an industry-by-industry or industry-quarter by industry-quarter basis.  To save some computation time, I'll run it all in one go.  This is fine because the treatment propensity simulation above does not use industry-specific parameters to calculate treatment likelihoods."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "4d5cc0e3",
   "metadata": {
    "kernel": "Stata"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "\n",
      "Iteration 0:   log likelihood = -11341.123  \n",
      "Iteration 1:   log likelihood = -8458.2127  \n",
      "Iteration 2:   log likelihood = -7160.0203  \n",
      "Iteration 3:   log likelihood = -7088.6004  \n",
      "Iteration 4:   log likelihood = -7088.2951  \n",
      "Iteration 5:   log likelihood = -7088.2951  \n",
      "\n",
      "Logistic regression                                    Number of obs =  40,497\n",
      "                                                       LR chi2(1)    = 8505.66\n",
      "                                                       Prob > chi2   =  0.0000\n",
      "Log likelihood = -7088.2951                            Pseudo R2     =  0.3750\n",
      "\n",
      "------------------------------------------------------------------------------\n",
      "         tau | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]\n",
      "-------------+----------------------------------------------------------------\n",
      "           x |   -1.99339   .0298281   -66.83   0.000    -2.051852   -1.934928\n",
      "       _cons |  -.0128757   .0311895    -0.41   0.680    -.0740061    .0482546\n",
      "------------------------------------------------------------------------------\n",
      "\n",
      "(39,503 missing values generated)\n",
      "\n",
      "(39,503 missing values generated)\n"
     ]
    }
   ],
   "source": [
    "// logit\n",
    "logit tau x\n",
    "\n",
    "// store probabilities\n",
    "predict phat, pr\n",
    "\n",
    "// store regression weight\n",
    "gen wt = tau/phat + (1-tau)/(1-phat)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a9050db1",
   "metadata": {
    "kernel": "Stata"
   },
   "source": [
    "Impose overlap with a maxima-minima rule for estimated treatment probabilities."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "0271eb4a",
   "metadata": {
    "kernel": "Stata"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "\n",
      "\n",
      "\n",
      "\n",
      "\n",
      "\n",
      "(4,770 real changes made, 4,770 to missing)\n"
     ]
    }
   ],
   "source": [
    "// overlap\n",
    "qui su phat if tau == 0\n",
    "local min0 = `r(min)'\n",
    "local max0 = `r(max)'\n",
    "qui su phat if tau == 1\n",
    "local min1 = `r(min)'\n",
    "local max1 = `r(max)'\n",
    "replace wt = . if phat < max(`min0',`min1') | phat > min(`max0',`max1')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6b31aa42",
   "metadata": {
    "kernel": "Stata"
   },
   "source": [
    "Now, *within* a particular industry-quarter, the parameter $\\rho$ is fixed.  This allows us to freely estimate $E[y(\\tau=1,\\rho)|\\rho]$ and $E[\\tau=0,\\rho)|\\rho]$.  We will include $x$ as a covariate in the industry-quarter by industry-quarter regressions so that the expected outcomes we estimate are IPWRA estimators and thus doubly-robust."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "c47bbc07",
   "metadata": {
    "kernel": "Stata"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "\n",
      "\n",
      "\n",
      "\n",
      "\n",
      "(2,700 missing values generated)\n",
      "\n",
      "(2,700 missing values generated)\n"
     ]
    }
   ],
   "source": [
    "// within-group expected outcomes\n",
    "qui gen y1 = .\n",
    "qui gen y0 = .\n",
    "egen grp = group(j t)\n",
    "qui levelsof grp, local(glevels)\n",
    "foreach g of local glevels{\n",
    "    qui{\n",
    "    cap reg y tau [pw = wt] if grp == `g'\n",
    "    if _rc ~= 0{\n",
    "        continue\n",
    "    }\n",
    "    gen tau_orig = tau\n",
    "    replace tau = 1\n",
    "    predict y1tmp if e(sample)\n",
    "    replace tau = 0\n",
    "    predict y0tmp if e(sample)\n",
    "    replace tau = tau_orig\n",
    "    replace y1 = y1tmp if grp == `g'\n",
    "    replace y0 = y0tmp if grp == `g'\n",
    "    drop tau_orig y*tmp\n",
    "    }\n",
    "}\n",
    "// mean group outcomes\n",
    "bys grp: egen Ey0 = mean(y0)\n",
    "bys grp: egen Ey1 = mean(y1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f6480e78",
   "metadata": {
    "kernel": "Stata"
   },
   "source": [
    "To estimate the dose response functions $E[y(1,\\rho)]$ and $E[y(0,\\rho)]$ over industry-quarter level data, one could either use $\\rho$ as the continous treatment variable or a standardized form of $\\rho$.  In the case of the former, because it is bounded over $[0,1]$ and likely to be skewed, special care needs to be taken with specifying the link function used in the dose response estimation.  It turns out that it is more effective to standardize $\\rho$ (done below) and then un-standardize it post-estimation because of the good performance that we get out of the usual Gaussian distribution with identity link function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "b6807ded",
   "metadata": {
    "kernel": "Stata"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "\n",
      "                             rho\n",
      "-------------------------------------------------------------\n",
      "      Percentiles      Smallest\n",
      " 1%          .04              0\n",
      " 5%           .1              0\n",
      "10%          .16              0       Obs               1,600\n",
      "25%           .3              0       Sum of wgt.       1,600\n",
      "\n",
      "50%          .54                      Mean             .53455\n",
      "                        Largest       Std. dev.      .2712384\n",
      "75%          .76            .98\n",
      "90%           .9            .98       Variance       .0735703\n",
      "95%          .96              1       Skewness      -.0927423\n",
      "99%          .98              1       Kurtosis       1.900762\n",
      "\n",
      "\n",
      "\n",
      "       Transform |          k    [95% conf. interval]       Skewness\n",
      "-----------------+--------------------------------------------------\n",
      "      ln(-rho-k) |  -4.434525      (not calculated)        -.0000362\n",
      "\n",
      "\n",
      "\n",
      "\n",
      "\n",
      "(80,000 real changes made)\n",
      "\n",
      "\n",
      "                      ln(-rho+4.434525)\n",
      "-------------------------------------------------------------\n",
      "      Percentiles      Smallest\n",
      " 1%    -1.709036      -1.792517\n",
      " 5%    -1.626036      -1.792517\n",
      "10%    -1.379873      -1.709036       Obs               1,600\n",
      "25%    -.8213747      -1.709036       Sum of wgt.       1,600\n",
      "\n",
      "50%     .0146534                      Mean           4.94e-09\n",
      "                        Largest       Std. dev.      1.000306\n",
      "75%     .8744471       1.881574\n",
      "90%     1.353232       1.881574       Variance       1.000613\n",
      "95%     1.553642       1.881574       Skewness      -.0000364\n",
      "99%     1.751296       1.881574       Kurtosis       1.888488\n"
     ]
    }
   ],
   "source": [
    "// normalize rho\n",
    "su rho if _first, d\n",
    "\n",
    "cap drop lrho\n",
    "lnskew0 rho_deskewed = rho\n",
    "scalar lnskew0_k = `r(gamma)'\n",
    "qui su rho_deskewed, d\n",
    "scalar lnskew_sd = `r(sd)'\n",
    "scalar lnskew_mu = `r(mean)'\n",
    "replace rho_deskewed = (rho_deskewed-lnskew_mu)/lnskew_sd\n",
    "\n",
    "su rho_deskewed if _first, d"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4eea14ab",
   "metadata": {
    "kernel": "Stata"
   },
   "source": [
    "We can now estimate a dose response function on $E[y(0,\\rho)]$ with covariate $w$.  Recall that $x$ predicts firm-level treatment and was correlated with $w$, so $w$ should predict industry-quarter level treatment $\\rho$.  To speed things along, I've specified `flag(0)` to turn off covariate balancing tests (since we know that, by construction, `w` is sufficient to balance over $\\rho$)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "2aafd123",
   "metadata": {
    "kernel": "Stata",
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "file C:/Users/james/.stata_kernel_cache/graph0.svg saved as SVG format\n",
      "file C:/Users/james/.stata_kernel_cache/graph0.pdf saved as PDF format\n"
     ]
    }
   ],
   "source": [
    "// drf on y0\n",
    "qui{\n",
    "cap restore\n",
    "preserve\n",
    "keep if _first == 1\n",
    "drf w, outcome(Ey0) treatment(rho_deskewed) method(iwkernel) npoints(100) common(1) flag(0)\n",
    "matrix x = e(tvec)\n",
    "local npoints `=rowsof(x)'\n",
    "matrix y0_hat = e(b)\n",
    "mat y0_hat = y0_hat[1,1..`npoints']'\n",
    "svmat x, name(rho_points)\n",
    "svmat y0_hat, name(y0_points)\n",
    "replace rho_points1 = -1*lnskew0_k-exp(rho_points1*lnskew_sd+lnskew_mu)\n",
    "gen ytrue = b0 + b3*log(.1+rho_points1)\n",
    "twoway (scatter y0_points rho_points1) (scatter ytrue rho_points1)\n",
    "restore\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "76268b80",
   "metadata": {
    "kernel": "Stata"
   },
   "source": [
    "![stata plot](https://nordlund.ai/notebooks/assets/sutnva0.svg)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "81a0d21a",
   "metadata": {
    "kernel": "Stata"
   },
   "source": [
    "We do the same for $E[y(1,\\rho)]$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "06152e5c",
   "metadata": {
    "kernel": "Stata"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "file C:/Users/james/.stata_kernel_cache/graph1.svg saved as SVG format\n",
      "file C:/Users/james/.stata_kernel_cache/graph1.pdf saved as PDF format\n"
     ]
    }
   ],
   "source": [
    "// drf on y1\n",
    "qui{\n",
    "preserve\n",
    "keep if _first == 1\n",
    "drf w, outcome(Ey1) treatment(rho_deskewed) method(iwkernel) npoints(100) common(1) flag(0)\n",
    "matrix x = e(tvec)\n",
    "local npoints `=rowsof(x)'\n",
    "matrix y1_hat = e(b)\n",
    "mat y1_hat = y1_hat[1,1..`npoints']'\n",
    "svmat x, name(rho_points)\n",
    "svmat y1_hat, name(y1_points)\n",
    "replace rho_points1 = -1*lnskew0_k-exp(rho_points1*lnskew_sd+lnskew_mu)\n",
    "gen ytrue = b0 + b1 + b2*log(.1+rho_points1)\n",
    "twoway (scatter y1_points rho_points1) (scatter ytrue rho_points1)\n",
    "restore\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "88a0a3e3",
   "metadata": {
    "kernel": "Stata"
   },
   "source": [
    "![stata graph](https://nordlund.ai/notebooks/assets/sutnva1.svg)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "49076147",
   "metadata": {
    "kernel": "Stata"
   },
   "source": [
    "The estimated effects for $\\tau=0$ and for $\\tau=1$ firms line up very well with their true values.  The only failure is at the tail end of the distribution for $\\rho$.  Given that the dose response function is estimated nonparameterically, failures near the edge of the distribution are expected."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Stata",
   "language": "stata",
   "name": "stata"
  },
  "language_info": {
   "codemirror_mode": "stata",
   "file_extension": ".do",
   "mimetype": "text/x-stata",
   "name": "stata",
   "version": "15.1"
  },
  "sos": {
   "kernels": [
    [
     "Stata",
     "stata",
     "Stata",
     "#CAE8F3",
     ""
    ]
   ],
   "panel": {
    "displayed": true,
    "height": 0
   },
   "version": "0.22.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}