*! causalimpact 1.0.0  05aug2026
*! Bayesian structural time-series causal impact analysis
*! Author: Dr Merwan Roudane  <merwanroudane920@gmail.com>
*! GitHub: https://github.com/merwanroudane
*!
*! Implements Brodersen, Gallusser, Koehler, Remy and Scott (2015),
*! "Inferring causal impact using Bayesian structural time-series models",
*! Annals of Applied Statistics 9(1), 247-274.  DOI: 10.1214/14-AOAS788
*! Compatible with the R package CausalImpact 1.4.1 (google/CausalImpact),
*! whose source is bundled under reference_R/ for step-by-step comparison.
*!
*! ---------------------------------------------------------------------------
*! STEP -> EQUATION MAP        (full derivation: help causalimpact_methods)
*! ---------------------------------------------------------------------------
*!  S1  input / period checks ............ Sec. 1; R FormatInputForCausalImpact
*!  S2  pre-period standardisation ....... R impact_misc.R Standardize()
*!  S3  post-period response set missing . Sec. 2.3; R impact_analysis.R:411
*!  S4  state-space assembly ............. eq. (2.1)-(2.2)
*!  S5  local level block ................ eq. (2.3) with delta_t == 0
*!  S6  seasonal block ................... eq. (2.5)
*!  S7  static regression block .......... Sec. 2.1 "static coefficients"
*!  S8  dynamic regression block ......... eq. (2.6)
*!  S9  variance priors .................. eq. (2.7)
*! S10  spike-and-slab prior ............. eq. (2.8)-(2.9)
*! S11  slab / Zellner g-prior ........... eq. (2.10)-(2.12)
*! S12  Gibbs state draw (Durbin-Koopman)  Sec. 2.3 "Posterior simulation"
*! S13  Gibbs gamma / beta / sigma draw .. eq. (2.13)
*! S14  posterior predictive counterfactual eq. (2.14)
*! S15  pointwise causal effect .......... eq. (2.15)
*! S16  cumulative causal effect ......... eq. (2.16)
*! S17  running-average causal effect .... eq. (2.17)
*! S18  summary table / tail-area p ...... R CompileSummaryTable()
*! S19  verbal report .................... R InterpretSummaryTable()
*! S20  three-panel plot ................. Fig. 1, Fig. 5-7; R impact_plot.R
*! ---------------------------------------------------------------------------

program define causalimpact, eclass sortpreserve
    version 14.0

    * Replay.  The options that control DISPLAY must be honoured here too,
    * otherwise `causalimpact, report' silently redisplays the table only.
    if replay() {
        if ("`e(cmd)'" != "causalimpact") error 301
        syntax [, REPort NOTABle DIgits(integer 2) ]
        if ("`notable'" == "") _ci_display
        if ("`report'"  != "") _ci_report, digits(`digits')
        exit
    }

    syntax varlist(numeric ts min=1) [if] [in] ,        ///
        PREperiod(numlist min=2 max=2)                  ///
        POSTperiod(numlist min=2 max=2)                 ///
        [                                               ///
          NITer(integer 1000)                           ///
          BURNin(real -1)                               ///
          SEED(string)                                  ///
          ALpha(real 0.05)                              ///
          Level(cilevel)                                ///
          PRIORLevelsd(real 0.01)                       ///
          NSEasons(integer 1)                           ///
          SEASONDuration(integer 1)                     ///
          DYNamicregression                             ///
          NOSTANDardize                                 ///
          MAXFlips(integer -1)                          ///
          MODELsize(real 3)                             ///
          R2(real 0.8)                                  ///
          PRIORDf(real 50)                              ///
          GINFo(real 0.01)                              ///
          DSHRinkage(real 0.5)                          ///
          SIGMAUpper(real -1)                           ///
          NOCONstant                                    ///
          GRaph                                         ///
          METRics(string)                               ///
          NAME(string)                                  ///
          SAVing(string)                                ///
          LEGend                                        ///
          COEFplot                                      ///
          REPort                                        ///
          DIgits(integer 2)                             ///
          GENerate(string)                              ///
          REPLACE                                       ///
          NOTABle                                       ///
          NODOTs                                        ///
        ]

    * ---------------------------------------------------------------- S1 ---
    qui tsset
    local tvar   "`r(timevar)'"
    local pvar   "`r(panelvar)'"
    local tdelta = r(tdelta)
    if ("`tvar'" == "") {
        di as err "data must be {help tsset:tsset} before using causalimpact"
        exit 459
    }
    if ("`pvar'" != "") {
        di as err "causalimpact analyses a single time series, but the data are"
        di as err "{help xtset:xtset} on panel variable {bf:`pvar'}."
        di as err "Select one unit, e.g.  causalimpact y x1 if `pvar'==1, ..."
        exit 459
    }

    if (`level' != c(level)) local alpha = 1 - `level'/100
    else                     local level = 100*(1 - `alpha')
    if (`alpha' <= 0 | `alpha' >= 1) {
        di as err "alpha() must lie strictly between 0 and 1"
        exit 198
    }

    tsrevar `varlist'
    local rvlist "`r(varlist)'"
    gettoken depvar  covars  : rvlist
    gettoken dispdep dispcov : varlist

    local J : word count `covars'
    local hasx = (`J' > 0)

    marksample touse, novarlist
    if (`hasx') markout `touse' `covars'
    qui count if `touse'
    if (r(N) < 4) {
        di as err "too few usable observations (need at least 4)"
        exit 2001
    }

    * ---------------------------------------------------------------- S1 ---
    if (`niter' < 10) {
        di as err "must draw, at the very least, 10 MCMC samples; recommending 1000"
        exit 198
    }
    if (`niter' < 1000) {
        di as txt "{p 0 2}note: results potentially inaccurate; consider using" ///
                  " more MCMC samples, e.g. niter(1000).{p_end}"
    }
    if (`burnin' < 0) local burnin = round(0.1*`niter')
    local burnin = round(`burnin')
    if (`burnin' < 1) local burnin = 1
    if (`burnin' >= `niter') {
        di as err "burnin() must be smaller than niter()"
        exit 198
    }
    if (`priorlevelsd' <= 0) {
        di as err "priorlevelsd() must be strictly positive"
        exit 198
    }
    if (`nseasons' < 1) {
        di as err "nseasons() cannot be 0; use 1 in order not to have a seasonal component"
        exit 198
    }
    if (`seasonduration' < 1) {
        di as err "seasonduration() must be at least 1"
        exit 198
    }
    if (`maxflips' <= 0 & `maxflips' != -1) {
        di as err "maxflips() must be positive, or -1 for no limit"
        exit 198
    }
    if (`r2' <= 0 | `r2' >= 1) {
        di as err "r2() must lie strictly between 0 and 1"
        exit 198
    }
    if (`priordf' <= 0) {
        di as err "priordf() must be strictly positive"
        exit 198
    }
    if (`modelsize' <= 0) {
        di as err "modelsize() must be strictly positive"
        exit 198
    }
    if (`dshrinkage' < 0 | `dshrinkage' > 1) {
        di as err "dshrinkage() must lie in [0,1]"
        exit 198
    }
    if (`ginfo' <= 0) {
        di as err "ginfo() must be strictly positive"
        exit 198
    }

    local dynreg = ("`dynamicregression'" != "")
    local stdz   = ("`nostandardize'" == "")
    local addcon = ("`noconstant'" == "")
    local dots   = ("`nodots'" == "")
    if (`dynreg' & !`hasx') {
        di as err "dynamicregression requires at least one covariate"
        exit 198
    }

    if ("`seed'" != "") {
        capture set seed `seed'
        if (_rc) {
            di as err "invalid seed()"
            exit 198
        }
    }

    gettoken pre1  pre2  : preperiod
    gettoken post1 post2 : postperiod
    if (`pre1' > `pre2') {
        di as err "preperiod(): the first element must not exceed the second"
        exit 198
    }
    if (`post1' > `post2') {
        di as err "postperiod(): the first element must not exceed the second"
        exit 198
    }
    if (`post1' <= `pre2') {
        di as err "the post-period must start after the pre-period ends"
        exit 198
    }

    sort `tvar'

    tempvar gapchk
    qui gen double `gapchk' = `tvar' - `tvar'[_n-1] if `touse' & `touse'[_n-1]
    qui summarize `gapchk' if `touse', meanonly
    local maxgap = r(max)
    if (`maxgap' < . & `tdelta' < . & `maxgap' > `tdelta') {
        di as txt "{p 0 2}note: the selected sample has gaps in {bf:`tvar'}." ///
                  " Like the R package, causalimpact models the data in row"  ///
                  " order and ignores the gaps; consider {help tsfill} if a"  ///
                  " seasonal component is used.{p_end}"
    }

    * ------------------------------------------------------- S2 .. S17 -----
    tempvar s_resp s_cumresp s_pred s_predlo s_predup s_cpred s_cpredlo
    tempvar s_cpredup s_eff s_efflo s_effup s_ceff s_cefflo s_ceffup s_avgeff
    foreach v in s_resp s_cumresp s_pred s_predlo s_predup s_cpred s_cpredlo ///
                 s_cpredup s_eff s_efflo s_effup s_ceff s_cefflo s_ceffup    ///
                 s_avgeff {
        qui gen double ``v'' = .
    }

    tempname SUM INCL BB VV

    local e_ci_ok  = 1
    local e_ci_msg ""

    mata: _ci_main()

    if (`e_ci_ok' == 0) {
        di as err "`e_ci_msg'"
        di as err "inference aborted"
        exit 498
    }

    * ---------------------------------------------------------------- S20 ---
    if ("`generate'" != "") {
        local gvars pred lower upper cum_pred cum_lower cum_upper           ///
                    effect effect_lower effect_upper cum_effect             ///
                    cum_effect_lower cum_effect_upper response cum_response ///
                    avg_effect
        local gsrc  `s_pred' `s_predlo' `s_predup' `s_cpred' `s_cpredlo'    ///
                    `s_cpredup' `s_eff' `s_efflo' `s_effup' `s_ceff'        ///
                    `s_cefflo' `s_ceffup' `s_resp' `s_cumresp' `s_avgeff'
        if ("`replace'" != "") {
            foreach g of local gvars {
                capture drop `generate'_`g'
            }
        }
        foreach g of local gvars {
            capture confirm new variable `generate'_`g'
            if (_rc) {
                di as err "variable `generate'_`g' already exists; use the replace option"
                exit 110
            }
        }
        local k = 1
        foreach g of local gvars {
            local src : word `k' of `gsrc'
            qui gen double `generate'_`g' = `src'
            local k = `k' + 1
        }
        label variable `generate'_pred       "Counterfactual (posterior mean)"
        label variable `generate'_lower      "Counterfactual lower `level'% bound"
        label variable `generate'_upper      "Counterfactual upper `level'% bound"
        label variable `generate'_effect     "Pointwise causal effect"
        label variable `generate'_cum_effect "Cumulative causal effect"
        label variable `generate'_avg_effect "Running-average causal effect"
    }

    * Graphs are built BEFORE any matrix is moved into e()  (stata-package #6)
    if ("`graph'" != "" | "`saving'" != "" | "`metrics'" != "") {
        _ci_plot, tvar(`tvar') touse(`touse')                                   ///
            resp(`s_resp') pred(`s_pred') predlo(`s_predlo') predup(`s_predup') ///
            eff(`s_eff') efflo(`s_efflo') effup(`s_effup')                      ///
            ceff(`s_ceff') cefflo(`s_cefflo') ceffup(`s_ceffup')                ///
            pre2(`pre2') post2(`post2')                                         ///
            metrics(`metrics') name(`name') saving(`saving')                    ///
            level(`level') depvar(`dispdep') `legend'
    }
    if ("`coefplot'" != "" & `hasx') {
        _ci_coefplot, incl(`INCL') name(`name')
    }

    * ------------------------------------------------------------ ereturn ---
    if (`hasx') {
        ereturn post `BB' `VV', esample(`touse')
    }
    else {
        ereturn post, esample(`touse')
    }
    ereturn local depvar "`dispdep'"

    ereturn matrix summary = `SUM'
    if (`hasx') ereturn matrix inclusion = `INCL'

    ereturn scalar N                = `e_ci_N'
    ereturn scalar N_pre            = `e_ci_npre'
    ereturn scalar N_post           = `e_ci_npost'
    ereturn scalar N_model          = `e_ci_nmodel'
    ereturn scalar niter            = `niter'
    ereturn scalar burnin           = `burnin'
    ereturn scalar ndraws           = `niter' - `burnin'
    ereturn scalar alpha            = `alpha'
    ereturn scalar level            = `level'
    ereturn scalar p                = `e_ci_p'
    ereturn scalar prob_effect      = 1 - `e_ci_p'
    ereturn scalar actual_avg       = `e_ci_act1'
    ereturn scalar actual_cum       = `e_ci_act2'
    ereturn scalar pred_avg         = `e_ci_pred1'
    ereturn scalar pred_cum         = `e_ci_pred2'
    ereturn scalar abseffect        = `e_ci_abs1'
    ereturn scalar abseffect_sd     = `e_ci_abssd1'
    ereturn scalar abseffect_cum    = `e_ci_abs2'
    ereturn scalar abseffect_cum_sd = `e_ci_abssd2'
    ereturn scalar releffect        = `e_ci_rel1'
    ereturn scalar releffect_sd     = `e_ci_relsd1'
    ereturn scalar sigma_obs        = `e_ci_sobs'
    ereturn scalar sigma_level      = `e_ci_slev'
    if (`nseasons' > 1) ereturn scalar sigma_seas = `e_ci_sseas'
    ereturn scalar priorlevelsd     = `priorlevelsd'
    ereturn scalar nseasons         = `nseasons'
    ereturn scalar seasonduration   = `seasonduration'
    ereturn scalar modelsize        = `modelsize'
    ereturn scalar r2               = `r2'
    ereturn scalar priordf          = `priordf'
    ereturn scalar ginfo            = `ginfo'
    ereturn scalar dshrinkage       = `dshrinkage'
    ereturn scalar maxflips         = `maxflips'
    ereturn scalar k_covariates     = `J'

    local regtxt "none"
    if (`hasx' & `dynreg')  local regtxt "dynamic (time-varying coefficients)"
    if (`hasx' & !`dynreg') local regtxt "static spike-and-slab"
    local modtxt "local level"
    if (`nseasons' > 1) local modtxt "`modtxt' + seasonal(`nseasons')"
    if (`hasx')         local modtxt "`modtxt' + regression"
    local stdtxt "no"
    if (`stdz') local stdtxt "yes"

    ereturn local pre_period  "`pre1' `pre2'"
    ereturn local post_period "`post1' `post2'"
    ereturn local timevar     "`tvar'"
    ereturn local covariates  "`dispcov'"
    ereturn local standardize "`stdtxt'"
    ereturn local regression  "`regtxt'"
    ereturn local model       "`modtxt'"
    ereturn local predict     "causalimpact_p"
    ereturn local title       "Bayesian structural time-series causal impact"
    ereturn local cmdline     "causalimpact `0'"
    ereturn local cmd         "causalimpact"

    if ("`notable'" == "") _ci_display
    if ("`report'"  != "") _ci_report, digits(`digits')
end


* =========================================================================== *
*  S18  Summary table                                                          *
* =========================================================================== *
program define _ci_display
    version 14.0

    tempname S
    matrix `S' = e(summary)

    local lev = e(level)
    local ci  = "`lev'% CI"
    local p   = e(p)

    local act1 = `S'[1,1]
    local act2 = `S'[2,1]
    local pr1  = `S'[1,2]
    local pr2  = `S'[2,2]
    local pl1  = `S'[1,3]
    local pl2  = `S'[2,3]
    local pu1  = `S'[1,4]
    local pu2  = `S'[2,4]
    local ps1  = `S'[1,5]
    local ps2  = `S'[2,5]
    local ae1  = `S'[1,6]
    local ae2  = `S'[2,6]
    local al1  = `S'[1,7]
    local al2  = `S'[2,7]
    local au1  = `S'[1,8]
    local au2  = `S'[2,8]
    local as1  = `S'[1,9]
    local as2  = `S'[2,9]
    local re1  = 100*`S'[1,10]
    local re2  = 100*`S'[2,10]
    local rl1  = 100*`S'[1,11]
    local rl2  = 100*`S'[2,11]
    local ru1  = 100*`S'[1,12]
    local ru2  = 100*`S'[2,12]
    local rs1  = 100*`S'[1,13]
    local rs2  = 100*`S'[2,13]

    local star " "
    if (`rl1' > 0 | `ru1' < 0) local star "*"
    if (`p' < 0.05)            local star "**"
    if (`p' < 0.01)            local star "***"

    local cov = trim("`e(covariates)'")
    if ("`cov'" == "") local cov "(none)"
    if (length("`cov'") > 29) local cov = substr("`cov'", 1, 26) + "..."
    local mdl = trim("`e(model)'")
    if (length("`mdl'") > 29) local mdl = substr("`mdl'", 1, 26) + "..."
    local reg = trim("`e(regression)'")
    if (length("`reg'") > 29) local reg = substr("`reg'", 1, 26) + "..."
    local dep = abbrev("`e(depvar)'", 29)

    di ""
    di as txt "{hline 79}"
    di as txt "Bayesian causal impact analysis  (Brodersen et al. 2015, AOAS 9:247-274)"
    di as txt "{hline 79}"
    di as txt "Response         " as res "`dep'"                             ///
       as txt _col(48) "Obs (total)  =" as res %10.0f e(N)
    di as txt "Covariates       " as res "`cov'"                             ///
       as txt _col(48) "Obs (pre)    =" as res %10.0f e(N_pre)
    di as txt "State model      " as res "`mdl'"                             ///
       as txt _col(48) "Obs (post)   =" as res %10.0f e(N_post)
    di as txt "Regression       " as res "`reg'"                             ///
       as txt _col(48) "Covariates   =" as res %10.0f e(k_covariates)
    di as txt "Standardised     " as res "`e(standardize)'"                  ///
       as txt _col(48) "MCMC draws   =" as res %10.0f e(niter)
    di as txt "Pre-period       " as res "`e(pre_period)'"                   ///
       as txt _col(48) "Burn-in      =" as res %10.0f e(burnin)
    di as txt "Post-period      " as res "`e(post_period)'"                  ///
       as txt _col(48) "Retained     =" as res %10.0f e(ndraws)
    di as txt "{hline 79}"
    di as txt _col(32) "Average" _col(58) "Cumulative"
    di as txt "{hline 79}"
    di as txt "Actual"                as res _col(28) %12.4g `act1' _col(53) %14.4g `act2'
    di as txt "Prediction"            as res _col(28) %12.4g `pr1'  _col(53) %14.4g `pr2'
    di as txt "  posterior s.d."      as res _col(28) %12.4g `ps1'  _col(53) %14.4g `ps2'
    di as txt "  `ci'"                as res _col(24) "[" %9.4g `pl1' "," %9.4g `pu1' "]" ///
                                             _col(50) "[" %10.4g `pl2' "," %10.4g `pu2' "]"
    di as txt "{hline 79}"
    di as txt "Absolute effect"       as res _col(28) %12.4g `ae1'  _col(53) %14.4g `ae2'
    di as txt "  posterior s.d."      as res _col(28) %12.4g `as1'  _col(53) %14.4g `as2'
    di as txt "  `ci'"                as res _col(24) "[" %9.4g `al1' "," %9.4g `au1' "]" ///
                                             _col(50) "[" %10.4g `al2' "," %10.4g `au2' "]"
    di as txt "{hline 79}"
    di as txt "Relative effect (%)"   as res _col(28) %12.3f `re1'  _col(53) %14.3f `re2'
    di as txt "  posterior s.d."      as res _col(28) %12.3f `rs1'  _col(53) %14.3f `rs2'
    di as txt "  `ci'"                as res _col(24) "[" %9.3f `rl1' "," %9.3f `ru1' "]" ///
                                             _col(50) "[" %10.3f `rl2' "," %10.3f `ru2' "]"
    di as txt "{hline 79}"
    di as txt "Posterior tail-area probability  p" as res _col(40) %12.5f `p' "  " as txt "`star'"
    di as txt "Posterior probability of an effect" as res _col(40) %11.4f 100*(1-`p') "%"
    di as txt "{hline 79}"
    di as txt "*** p<.01, ** p<.05, * `lev'% credible interval of the effect excludes 0"
    di as txt "{p 0 2}Intervals are central (1-alpha) posterior credible intervals of" ///
              " the counterfactual; p is the one-sided Bayesian tail-area probability" ///
              " of the cumulative effect (eq. 2.16).{p_end}"

    if (e(k_covariates) > 0) _ci_incltable

    di as txt "{p 0 2}Type {stata causalimpact, report:causalimpact, report}" ///
              " for a verbal interpretation, or" ///
              " {stata help causalimpact_methods:help causalimpact methods}" ///
              " for the model.{p_end}"
    di ""
end


* =========================================================================== *
*  S10-S11  Spike-and-slab inclusion table                                     *
* =========================================================================== *
program define _ci_incltable
    version 14.0
    tempname I
    capture matrix `I' = e(inclusion)
    if (_rc) exit
    local Jd = rowsof(`I')
    local rn : rownames `I'

    di ""
    di as txt "Regression component: posterior summary of the synthetic control"
    di as txt "{hline 79}"
    di as txt "Covariate" _col(24) "P(incl)" _col(38) "Post. mean" ///
              _col(54) "Post. s.d." _col(68) "P(b>0)"
    di as txt "{hline 79}"
    forvalues j = 1/`Jd' {
        local nm : word `j' of `rn'
        di as txt %-22s abbrev("`nm'",22) as res             ///
           _col(23) %11.3f `I'[`j',1] _col(38) %12.4g `I'[`j',2] ///
           _col(54) %11.4g `I'[`j',3] _col(66) %11.3f `I'[`j',4]
    }
    di as txt "{hline 79}"
    di as txt "{p 0 2}P(incl) is the posterior inclusion probability of the" ///
              " spike-and-slab prior (eq. 2.8-2.9). Posterior mean and s.d." ///
              " are model-averaged: draws that exclude a covariate contribute" ///
              " a coefficient of exactly zero. P(b>0) is conditional on" ///
              " inclusion. As in the R package, coefficients are reported in" ///
              " the metric the model is fitted in (standardised, unless" ///
              " nostandardize was specified).{p_end}"
end


* =========================================================================== *
*  S19  Verbal report -- faithful port of R InterpretSummaryTable()            *
* =========================================================================== *
program define _ci_report
    version 14.0
    syntax [, DIgits(integer 2) ]

    if ("`e(cmd)'" != "causalimpact") error 301

    tempname S
    matrix `S' = e(summary)
    local alpha = e(alpha)
    local ci    = string(round(100*(1-`alpha'))) + "%"

    mata: _ci_report_strings(`digits')

    local sig = (!(`S'[1,11] < 0 & `S'[1,12] > 0))
    local pos = (`S'[1,10] > 0)
    local p   = e(p)
    * R prints round(p, 3) with a leading zero and trailing zeros dropped
    * ("0.001", "0.185"); Stata's string() drops the leading zero. Restore it
    * so the report is textually identical to summary(impact, "report") in R.
    local pr  = string(round(`p',0.001))
    if (substr("`pr'",1,1) == ".")  local pr = "0`pr'"
    if (substr("`pr'",1,2) == "-.") local pr = "-0" + substr("`pr'",2,.)

    local by1 "In"
    local by2 "Had"
    if (`sig') {
        local by1 "By contrast, in"
        local by2 "By contrast, had"
    }
    local dir "a decrease of"
    if (`pos') local dir "an increase of"

    di ""
    di as txt "{hline 79}"
    di as txt "Analysis report {causalimpact}"
    di as txt "{hline 79}"
    di ""

    _ci_wrap `"During the post-intervention period, the response variable had an average value of approx. `r_act1'. `by1' the absence of an intervention, we would have expected an average response of `r_pred1'. The `ci' interval of this counterfactual prediction is [`r_predlo1', `r_predup1']. Subtracting this prediction from the observed response yields an estimate of the causal effect the intervention had on the response variable. This effect is `r_abs1' with a `ci' interval of [`r_abslo1', `r_absup1']. For a discussion of the significance of this effect, see below."'
    di ""
    _ci_wrap `"Summing up the individual data points during the post-intervention period (which can only sometimes be meaningfully interpreted), the response variable had an overall value of `r_act2'. `by2' the intervention not taken place, we would have expected a sum of `r_pred2'. The `ci' interval of this prediction is [`r_predlo2', `r_predup2']."'
    di ""
    _ci_wrap `"The above results are given in terms of absolute numbers. In relative terms, the response variable showed `dir' `r_rel1'. The `ci' interval of this percentage is [`r_rello1', `r_relup1']."'

    local gap = abs(`S'[1,10] - `S'[1,6]/`S'[1,2])
    if (`gap' >= 0.05) {
        di ""
        _ci_wrap `"(Note that the expected relative effect isn't generally the same as the expected absolute effect divided by the expected prediction since distributions are often not symmetric.)"'
    }

    di ""
    if (`sig' & `pos') {
        _ci_wrap `"This means that the positive effect observed during the intervention period is statistically significant and unlikely to be due to random fluctuations. It should be noted, however, that the question of whether this increase also bears substantive significance can only be answered by comparing the absolute effect (`r_abs1') to the original goal of the underlying intervention."'
    }
    else if (`sig' & !`pos') {
        _ci_wrap `"This means that the negative effect observed during the intervention period is statistically significant. If the experimenter had expected a positive effect, it is recommended to double-check whether anomalies in the control variables may have caused an overly optimistic expectation of what should have happened in the response variable in the absence of the intervention."'
    }
    else if (!`sig' & `pos') {
        _ci_wrap `"This means that, although the intervention appears to have caused a positive effect, this effect is not statistically significant when considering the entire post-intervention period as a whole. Individual days or shorter stretches within the intervention period may of course still have had a significant effect, as indicated whenever the lower limit of the impact time series (lower plot) was above zero."'
    }
    else {
        _ci_wrap `"This means that, although it may look as though the intervention has exerted a negative effect on the response variable when considering the intervention period as a whole, this effect is not statistically significant, and so cannot be meaningfully interpreted."'
    }

    if (!`sig') {
        di ""
        _ci_wrap `"The apparent effect could be the result of random fluctuations that are unrelated to the intervention. This is often the case when the intervention period is very long and includes much of the time when the effect has already worn off. It can also be the case when the intervention period is too short to distinguish the signal from the noise. Finally, failing to find a significant effect can happen when there are not enough control variables or when these variables do not correlate well with the response variable during the learning period."'
    }

    di ""
    if (`p' < `alpha') {
        _ci_wrap `"The probability of obtaining this effect by chance is very small (Bayesian one-sided tail-area probability p = `pr'). This means the effect is statistically significant. It can be considered causal if the model assumptions are satisfied."'
    }
    else {
        _ci_wrap `"The probability of obtaining this effect by chance is p = `pr'. This means the effect may be spurious and would generally not be considered statistically significant."'
    }

    di ""
    _ci_wrap `"For more details, including the model assumptions behind the method, see help causalimpact_methods and Brodersen et al. (2015), Annals of Applied Statistics 9(1), 247-274."'
    di as txt "{hline 79}"
    di ""
end


program define _ci_wrap
    version 14.0
    args txt
    local w 78
    local s `"`txt'"'
    while (`"`s'"' != "") {
        local L = length(`"`s'"')
        if (`L' <= `w') {
            di as txt `"`s'"'
            local s ""
        }
        else {
            local cut = `w' + 1
            while (`cut' > 1 & substr(`"`s'"',`cut',1) != " ") {
                local cut = `cut' - 1
            }
            if (`cut' <= 1) local cut = `w' + 1
            di as txt `"`=substr(`"`s'"',1,`cut'-1)'"'
            local s = trim(substr(`"`s'"',`cut'+1,.))
        }
    }
end


* =========================================================================== *
*  S20  Three-panel plot -- Fig. 1 and Fig. 5-7 of the paper                    *
* =========================================================================== *
program define _ci_plot
    version 14.0
    syntax , tvar(varname) touse(varname)                                    ///
             resp(varname) pred(varname) predlo(varname) predup(varname)     ///
             eff(varname) efflo(varname) effup(varname)                      ///
             ceff(varname) cefflo(varname) ceffup(varname)                   ///
             pre2(string) post2(string)                                      ///
             level(string) depvar(string)                                    ///
             [ metrics(string) name(string) saving(string) LEGend ]

    if ("`metrics'" == "") local metrics "original pointwise cumulative"
    local metrics = lower("`metrics'")

    local gname "causalimpact"
    if ("`name'" != "") local gname "`name'"

    * ---- R palette ------------------------------------------------------
    * ggplot2 "slategray2"  = RGB 185 211 238   (the credible band)
    * ggplot2 "darkblue"    = RGB   0   0 139   (the counterfactual)
    * theme_bw panel border = grey20 = RGB  51  51  51
    * theme_bw grid lines   = grey92 = RGB 235 235 235
    * NOTE: RGB triples MUST be quoted in Stata, otherwise "185" is parsed
    * as a named style and silently falls back to the default colour.
    local cBAND  `""185 211 238""'
    local cPRED  `""0 0 139""'
    local cBORD  `""51 51 51""'
    local cGRID  `""235 235 235""'
    local cMARK  `""128 128 128""'

    local band  fcolor(`cBAND') lcolor(`cBAND') lwidth(none)
    local predl lcolor(`cPRED') lpattern(shortdash) lwidth(medthin)
    local respl lcolor(black) lwidth(medthin)

    * theme_bw look: white panel with a thin dark border and pale gridlines
    local sch   graphregion(color(white) margin(zero) lcolor(white))         ///
                plotregion(color(white) lcolor(`cBORD') lwidth(thin)         ///
                           margin(l=1 r=1 b=1 t=1))
    local ylab  ylabel(, labsize(vsmall) angle(0) format(%9.0g)              ///
                       grid glcolor(`cGRID') glwidth(vthin) glpattern(solid) ///
                       tlcolor(`cBORD'))
    local xlabS xlabel(, labsize(vsmall) nolabels                            ///
                       grid glcolor(`cGRID') glwidth(vthin) glpattern(solid) ///
                       tlcolor(`cBORD'))
    local xlabB xlabel(, labsize(vsmall)                                     ///
                       grid glcolor(`cGRID') glwidth(vthin) glpattern(solid) ///
                       tlcolor(`cBORD'))
    local vln   xline(`pre2' `post2', lpattern(dash) lcolor(`cMARK')         ///
                      lwidth(vthin))
    local zln   yline(0, lcolor(`cMARK') lwidth(vthin) lpattern(solid))

    * facet-strip label, right-hand side, as in facet_grid(metric ~ .)
    local strip1 subtitle("original", pos(3) orientation(rvertical)          ///
                    size(vsmall) box bcolor(`cGRID') bmargin(small)          ///
                    lcolor(`cBORD') lwidth(vthin))
    local strip2 subtitle("pointwise", pos(3) orientation(rvertical)         ///
                    size(vsmall) box bcolor(`cGRID') bmargin(small)          ///
                    lcolor(`cBORD') lwidth(vthin))
    local strip3 subtitle("cumulative", pos(3) orientation(rvertical)        ///
                    size(vsmall) box bcolor(`cGRID') bmargin(small)          ///
                    lcolor(`cBORD') lwidth(vthin))

    * R's plot carries no legend and no axis titles; both are opt-in here.
    local leg legend(off)
    if ("`legend'" != "") {
        local leg legend(order(3 "Observed" 2 "Counterfactual"               ///
                               1 "`level'% credible interval")               ///
                         rows(1) size(vsmall) region(lcolor(white))          ///
                         symxsize(6) pos(6))
    }

    local todraw ""
    local nsel = 0
    if (strpos("`metrics'","orig")) local nsel = `nsel' + 1
    if (strpos("`metrics'","point")) local nsel = `nsel' + 1
    if (strpos("`metrics'","cum")) local nsel = `nsel' + 1
    local kdone = 0

    if (strpos("`metrics'","orig")) {
        local kdone = `kdone' + 1
        local xl "`xlabS'"
        if (`kdone' == `nsel') local xl "`xlabB'"
        capture graph drop `gname'_original
        twoway (rarea `predup' `predlo' `tvar' if `touse', `band')       ///
               (line `pred' `tvar' if `touse', `predl')                  ///
               (line `resp' `tvar' if `touse', `respl'),                 ///
               `sch' `vln' `ylab' `xl' `strip1' `leg'                    ///
               ytitle("") xtitle("")                                     ///
               name(`gname'_original, replace) nodraw
        local todraw "`todraw' `gname'_original"
    }
    if (strpos("`metrics'","point")) {
        local kdone = `kdone' + 1
        local xl "`xlabS'"
        if (`kdone' == `nsel') local xl "`xlabB'"
        capture graph drop `gname'_pointwise
        twoway (rarea `effup' `efflo' `tvar' if `touse', `band')         ///
               (line `eff' `tvar' if `touse', `predl'),                  ///
               `sch' `vln' `zln' `ylab' `xl' `strip2'                    ///
               ytitle("") xtitle("") legend(off)                         ///
               name(`gname'_pointwise, replace) nodraw
        local todraw "`todraw' `gname'_pointwise"
    }
    if (strpos("`metrics'","cum")) {
        local kdone = `kdone' + 1
        capture graph drop `gname'_cumulative
        twoway (rarea `ceffup' `cefflo' `tvar' if `touse', `band')       ///
               (line `ceff' `tvar' if `touse', `predl'),                 ///
               `sch' `vln' `zln' `ylab' `xlabB' `strip3'                 ///
               ytitle("") xtitle("") legend(off)                         ///
               name(`gname'_cumulative, replace) nodraw
        local todraw "`todraw' `gname'_cumulative"
    }

    local nplot : word count `todraw'
    if (`nplot' == 0) exit

    if (`nplot' == 1) {
        local final : word 1 of `todraw'
        graph display `final'
    }
    else {
        capture graph drop `gname'
        graph combine `todraw', cols(1)                                  ///
            graphregion(color(white) margin(zero)) imargin(zero)         ///
            ysize(`=2.0*`nplot'') xsize(6.5)                             ///
            name(`gname', replace)
        local final "`gname'"
    }
    if (`"`saving'"' != "") {
        graph save `final' `saving', replace
    }
end


* =========================================================================== *
*  Inclusion plot (companion to R plot(bsts.model, "coefficients"))            *
* =========================================================================== *
program define _ci_coefplot
    version 14.0
    syntax , incl(name) [ name(string) ]

    * Reproduces BoomSpikeSlab::PlotMarginalInclusionProbabilities, which is
    * what plot(bsts.model, "coefficients") draws in R: covariates ordered by
    * posterior inclusion probability with the largest at the top, each bar
    * shaded by the conditional probability that the coefficient is positive
    * given that it is in the model.  R uses a white-to-black ramp; here the
    * same information is carried on a blue ramp running from a pale blue
    * (negative) to darkblue RGB 0 0 139 (positive) -- the top of the ramp is
    * exactly the colour of the counterfactual line in the impact figure, so
    * the two figures read as one family.  A single linear ramp is used rather
    * than one bent through slategray2, because the pale end and slategray2
    * are too close to tell apart at bar size.

    local gname "causalimpact_coef"
    if ("`name'" != "") local gname "`name'_coef"

    tempname I
    matrix `I' = `incl'
    local Jd = rowsof(`I')
    local rn : rownames `I'

    local cBORD `""51 51 51""'
    local cGRID `""235 235 235""'
    * ramp endpoints: pale blue -> darkblue (0 0 139), through slategray2
    local rLO = 222
    local gLO = 235
    local bLO = 250
    local rHI = 0
    local gHI = 0
    local bHI = 139

    local maxbar = 60
    local trunc  = 0
    if (`Jd' > `maxbar') {
        local trunc = 1
        di as txt "{p 0 2}note: only the `maxbar' covariates with the highest" ///
                  " posterior inclusion probability are plotted.{p_end}"
    }

    preserve
    clear
    qui set obs `Jd'
    qui svmat double `I', name(ci_)
    qui gen int ci_row = _n
    qui gen str32 ci_nm = ""
    forvalues j = 1/`Jd' {
        local nm : word `j' of `rn'
        qui replace ci_nm = "`nm'" in `j'
    }

    * largest inclusion probability at the TOP, as in R
    gsort ci_1 -ci_row
    qui gen int ci_ord = _n
    if (`trunc') {
        qui drop if ci_ord <= `Jd' - `maxbar'
        qui replace ci_ord = ci_ord - (`Jd' - `maxbar')
    }
    qui count
    local nb = r(N)

    local lb ""
    local cmd ""
    forvalues j = 1/`nb' {
        local nmj = ci_nm[`j']
        local lb  `"`lb' `j' "`=abbrev("`nmj'",18)'""'
        local pp  = ci_4[`j']
        if (`pp' >= .) local pp = 0.5
        if (`pp' < 0)  local pp = 0
        if (`pp' > 1)  local pp = 1
        local rr = round(`rLO' + `pp'*(`rHI' - `rLO'))
        local gg = round(`gLO' + `pp'*(`gHI' - `gLO'))
        local bb = round(`bLO' + `pp'*(`bHI' - `bLO'))
        local cbar `""`rr' `gg' `bb'""'
        local cmd `"`cmd' (bar ci_1 ci_ord if ci_ord==`j', horizontal"'
        local cmd `"`cmd' barwidth(0.68) fcolor(`cbar') lcolor(`cBORD')"'
        local cmd `"`cmd' lwidth(vthin))"'
    }

    capture graph drop `gname'
    twoway `cmd',                                                             ///
           graphregion(color(white) margin(zero))                             ///
           plotregion(color(white) lcolor(`cBORD') lwidth(thin))              ///
           title("Posterior inclusion probabilities", size(medsmall) pos(11)) ///
           subtitle("shading: P(coefficient > 0 | included)"                  ///
                    " -- pale negative, dark blue positive",                  ///
                    size(vsmall) pos(11) margin(b=2))                         ///
           xtitle("P(inclusion)", size(small))                                ///
           ytitle("")                                                         ///
           ylabel(`lb', angle(0) labsize(vsmall) nogrid                       ///
                  tlcolor(`cBORD'))                                           ///
           yscale(range(0.3 `=`nb'+0.7'))                                     ///
           xlabel(0(0.2)1, format(%3.1f) labsize(vsmall)                      ///
                  grid glcolor(`cGRID') glwidth(vthin) glpattern(solid))      ///
           xscale(range(0 1))                                                 ///
           legend(off)                                                        ///
           name(`gname', replace)
    restore
end


* =========================================================================== *
*                               MATA ENGINE                                    *
* =========================================================================== *
version 14.0
mata:

// Assembled state-space model, eq. (2.1)-(2.2)
struct cissm {
    real matrix    Tadv, Thold, P1
    real colvector a1, adv
    real scalar    d, S, dur, m, nx, dyn, dss
}

// -------------------------------------------------------------------------
// R-compatible quantile (stats::quantile type 7); xs sorted ascending.
// -------------------------------------------------------------------------
real scalar _ci_quantile(real colvector xs, real scalar p)
{
    real scalar n, h, lo
    n = rows(xs)
    if (n == 1) {
        return(xs[1])
    }
    h  = (n - 1)*p + 1
    lo = floor(h)
    if (lo < 1) {
        return(xs[1])
    }
    if (lo >= n) {
        return(xs[n])
    }
    return(xs[lo] + (h - lo)*(xs[lo+1] - xs[lo]))
}

real scalar _ci_qcol(real colvector x, real scalar p)
{
    real colvector xs
    xs = sort(x, 1)
    return(_ci_quantile(xs, p))
}

// log|A| via Cholesky; missing when A is not positive definite.
real scalar _ci_logdet(real matrix A)
{
    real matrix    L
    real colvector dg
    if (rows(A) == 0) {
        return(0)
    }
    L = cholesky(A)
    if (hasmissing(L)) {
        return(.)
    }
    dg = diagonal(L)
    if (min(dg) <= 0) {
        return(.)
    }
    return(2*sum(log(dg)))
}

// -------------------------------------------------------------------------
// S9  Draw sigma^2 from its inverse-Gamma full conditional, eq. (2.7),
//     truncated so that sigma <= uplim (Boom SdPrior upper.limit).
//       prior 1/sigma^2 ~ G(df/2, ss/2),  ss = df * guess^2
//       post  1/sigma^2 ~ G((df+n)/2, (ss+sumsq)/2)
// -------------------------------------------------------------------------
real scalar _ci_draw_sig2(real scalar df, real scalar ss, real scalar n,
                          real scalar sumsq, real scalar uplim)
{
    real scalar a, b, plow, u, prec, lo
    a = (df + n)/2
    b = 2/(ss + sumsq)
    if (a <= 0 | b <= 0) {
        return(.)
    }
    if (uplim >= . | uplim <= 0) {
        return(1/rgamma(1, 1, a, b))
    }
    lo   = 1/(uplim*uplim)
    plow = gammap(a, lo/b)
    if (plow > 0.9999999) {
        return(uplim*uplim)
    }
    u = plow + (1 - plow)*runiform(1, 1)
    if (u >= 1) {
        u = 1 - 1e-12
    }
    prec = b*invgammap(a, u)
    if (prec <= 0) {
        return(uplim*uplim)
    }
    return(1/prec)
}

// -------------------------------------------------------------------------
// S12  Kalman filter (univariate y, possibly missing).  Stores v, F, K,
//      which is all the fast state smoother needs.
// -------------------------------------------------------------------------
void _ci_filter(real colvector y, struct cissm scalar M, real matrix Zmat,
                real scalar s2obs, real colvector RQRd, real colvector a1use,
                real colvector v, real colvector F, real matrix K)
{
    real scalar    t, m, d, Ft
    real colvector a, Zt, PZ, Kt
    real matrix    P, Tt, L, RQ

    m  = M.m
    d  = M.d
    v  = J(m, 1, 0)
    F  = J(m, 1, .)
    K  = J(d, m, 0)
    a  = a1use
    P  = M.P1
    RQ = diag(RQRd)

    for (t = 1; t <= m; t++) {
        Tt = M.Thold
        if (M.adv[t] == 1) {
            Tt = M.Tadv
        }
        Zt = Zmat[., t]
        if (y[t] < .) {
            PZ = P*Zt
            Ft = Zt'PZ + s2obs
            if (Ft <= 0) {
                Ft = 1e-12
            }
            v[t]    = y[t] - Zt'a
            F[t]    = Ft
            Kt      = (Tt*PZ)/Ft
            K[., t] = Kt
            L       = Tt - Kt*Zt'
            a       = Tt*a + Kt*v[t]
            P       = Tt*P*L'
        }
        else {
            v[t] = 0
            F[t] = .
            a    = Tt*a
            P    = Tt*P*Tt'
        }
        P = P + RQ
        P = (P + P')/2
    }
}

// -------------------------------------------------------------------------
// S12  Fast state smoother of Durbin and Koopman (2002).
// -------------------------------------------------------------------------
real matrix _ci_smooth(struct cissm scalar M, real matrix Zmat,
                       real colvector v, real colvector F, real matrix K,
                       real colvector RQRd, real colvector a1use)
{
    real scalar    t, m, d
    real colvector r, Zt
    real matrix    ah, Tt, L, rstore

    m      = M.m
    d      = M.d
    ah     = J(d, m, 0)
    rstore = J(d, m + 1, 0)
    r      = J(d, 1, 0)

    for (t = m; t >= 1; t--) {
        Tt = M.Thold
        if (M.adv[t] == 1) {
            Tt = M.Tadv
        }
        Zt = Zmat[., t]
        if (F[t] < .) {
            L = Tt - K[., t]*Zt'
            r = Zt*(v[t]/F[t]) + L'r
        }
        else {
            r = Tt'r
        }
        rstore[., t] = r
    }

    ah[., 1] = a1use + M.P1*rstore[., 1]
    for (t = 1; t <= m - 1; t++) {
        Tt = M.Thold
        if (M.adv[t] == 1) {
            Tt = M.Tadv
        }
        ah[., t+1] = Tt*ah[., t] + RQRd:*rstore[., t+1]
    }
    return(ah)
}

// -------------------------------------------------------------------------
// S12  Durbin-Koopman (2002) simulation smoother.
//
// The textbook recipe alpha~ = ahat(y) - ahat(y+) + alpha+ needs TWO smoother
// passes.  Because the smoother is affine in (y, a1) with the same a1 in both,
//     ahat(y) - ahat(y+) = ahat(y - y+ ; a1 = 0),
// so a SINGLE pass on y* = y - y+ with a zero initial mean is algebraically
// identical at half the cost.  See help causalimpact_methods, Step 12.
// -------------------------------------------------------------------------
real matrix _ci_simsmooth(real colvector y, struct cissm scalar M,
                          real matrix Zmat, real scalar s2obs,
                          real colvector RQRd)
{
    real scalar    t, m, d
    real colvector ap, ystar, v, F, zero1, eta, sdq
    real matrix    aplus, K, ahat, Tt, Lch

    m     = M.m
    d     = M.d
    aplus = J(d, m, 0)
    ystar = J(m, 1, .)
    sdq   = sqrt(RQRd)

    Lch = cholesky(M.P1)
    if (hasmissing(Lch)) {
        Lch = diag(sqrt(diagonal(M.P1)))
    }
    ap = M.a1 + Lch*rnormal(d, 1, 0, 1)

    for (t = 1; t <= m; t++) {
        aplus[., t] = ap
        if (y[t] < .) {
            ystar[t] = y[t] - (Zmat[., t]'ap + rnormal(1, 1, 0, sqrt(s2obs)))
        }
        if (t < m) {
            Tt = M.Thold
            if (M.adv[t] == 1) {
                Tt = M.Tadv
            }
            eta = sdq:*rnormal(d, 1, 0, 1)
            ap  = Tt*ap + eta
        }
    }

    zero1 = J(d, 1, 0)
    _ci_filter(ystar, M, Zmat, s2obs, RQRd, zero1, v, F, K)
    ahat = _ci_smooth(M, Zmat, v, F, K, RQRd, zero1)

    return(ahat + aplus)
}

// -------------------------------------------------------------------------
// S13  Log marginal posterior of an inclusion vector, eq. (2.13):
//        V^-1 = X'X + Omega^-1
//        btil = V (X'ydot + Omega^-1 b)
//        N    = nu_eps + n
//        S    = s_eps + ydot'ydot + b'Omega^-1 b - btil' V^-1 btil
//      log p  = log p(gamma) + .5 log|Omega^-1| - .5 log|V^-1| - (N/2) log S
// -------------------------------------------------------------------------
real scalar _ci_logmarg(real colvector g, real matrix XtX, real colvector Xty,
                        real scalar yty, real matrix Om, real colvector bpr,
                        real scalar nu, real scalar seps, real scalar n,
                        real colvector logpi, real colvector log1mpi,
                        real scalar Sout)
{
    real colvector idx, xb, ob, btil
    real matrix    Omg, Vinv
    real scalar    k, ld1, ld2, Sg, Nn, lp

    idx = selectindex(g :> 0)
    k   = length(idx)
    Nn  = nu + n
    lp  = quadsum(g:*logpi) + quadsum((1 :- g):*log1mpi)

    if (k == 0) {
        Sg = seps + yty
        if (Sg <= 0) {
            return(.)
        }
        Sout = Sg
        return(lp - (Nn/2)*log(Sg))
    }
    if (k >= n) {
        return(.)
    }

    Omg  = Om[idx, idx]
    Vinv = XtX[idx, idx] + Omg
    ld1  = _ci_logdet(Omg)
    ld2  = _ci_logdet(Vinv)
    if (ld1 >= . | ld2 >= .) {
        return(.)
    }
    xb   = Xty[idx]
    ob   = Omg*bpr[idx]
    btil = cholsolve(Vinv, xb + ob)
    if (hasmissing(btil)) {
        btil = invsym(Vinv)*(xb + ob)
    }
    Sg = seps + yty + bpr[idx]'ob - btil'Vinv*btil
    if (Sg <= 0) {
        return(.)
    }
    Sout = Sg
    return(lp + 0.5*ld1 - 0.5*ld2 - (Nn/2)*log(Sg))
}

// -------------------------------------------------------------------------
// S13  One stochastic-search sweep over gamma, then conjugate draws of
//      sigma_eps^2 and beta.  g, beta and s2obs are updated in place.
// -------------------------------------------------------------------------
void _ci_ssvs(real colvector g, real colvector beta, real scalar s2obs,
              real matrix XtX, real colvector Xty, real scalar yty,
              real matrix Om, real colvector bpr, real scalar nu,
              real scalar seps, real scalar n, real colvector logpi,
              real colvector log1mpi, real colvector pivec,
              real scalar maxflips, real scalar sigupper)
{
    real scalar    Jd, j, jj, nflip, lp0, lp1, S0, S1, pr, Sg, Nn
    real colvector ordr, g0, g1, idx, xb, ob, btil, bs
    real matrix    Omg, Vinv, Vch

    Jd = rows(g)
    Nn = nu + n

    ordr = (1::Jd)
    ordr = ordr[order(runiform(Jd, 1), 1)]

    nflip = Jd
    if (maxflips > 0 & maxflips < Jd) {
        nflip = maxflips
    }

    for (jj = 1; jj <= nflip; jj++) {
        j = ordr[jj]
        if (pivec[j] >= 1) {
            g[j] = 1
            continue
        }
        if (pivec[j] <= 0) {
            g[j] = 0
            continue
        }
        g0 = g
        g1 = g
        g0[j] = 0
        g1[j] = 1
        S0 = .
        S1 = .
        lp0 = _ci_logmarg(g0, XtX, Xty, yty, Om, bpr, nu, seps, n, logpi, log1mpi, S0)
        lp1 = _ci_logmarg(g1, XtX, Xty, yty, Om, bpr, nu, seps, n, logpi, log1mpi, S1)
        if (lp0 >= . & lp1 >= .) {
            continue
        }
        if (lp1 >= .) {
            g[j] = 0
            continue
        }
        if (lp0 >= .) {
            g[j] = 1
            continue
        }
        if (lp0 > lp1) {
            pr = exp(lp1 - lp0)/(1 + exp(lp1 - lp0))
        }
        else {
            pr = 1/(1 + exp(lp0 - lp1))
        }
        g[j] = (runiform(1, 1) < pr)
    }

    Sg  = .
    lp0 = _ci_logmarg(g, XtX, Xty, yty, Om, bpr, nu, seps, n, logpi, log1mpi, Sg)
    if (Sg >= . | Sg <= 0) {
        Sg = seps + yty
    }

    s2obs = _ci_draw_sig2(Nn, Sg, 0, 0, sigupper)
    if (s2obs >= . | s2obs <= 0) {
        s2obs = Sg/Nn
    }

    beta = J(Jd, 1, 0)
    idx  = selectindex(g :> 0)
    if (length(idx) > 0) {
        Omg  = Om[idx, idx]
        Vinv = XtX[idx, idx] + Omg
        xb   = Xty[idx]
        ob   = Omg*bpr[idx]
        btil = cholsolve(Vinv, xb + ob)
        if (hasmissing(btil)) {
            btil = invsym(Vinv)*(xb + ob)
        }
        Vch = cholesky(invsym(Vinv))
        if (hasmissing(Vch)) {
            Vch = J(length(idx), length(idx), 0)
        }
        bs = btil + sqrt(s2obs)*(Vch*rnormal(length(idx), 1, 0, 1))
        beta[idx] = bs
    }
}

// -------------------------------------------------------------------------
// Ports of PrettifyNumber() / PrettifyPercentage() for the verbal report.
// -------------------------------------------------------------------------
string scalar _ci_prettynum(real scalar x, string scalar letter, real scalar rd)
{
    real scalar fz, ax
    if (x >= .) {
        return("NA")
    }
    ax = abs(x)
    if ((letter == "" & ax >= 1e9) | letter == "B") {
        return(strtrim(strofreal(x/1e9, "%30." + strofreal(rd) + "f")) + "B")
    }
    if ((letter == "" & ax >= 1e6) | letter == "M") {
        return(strtrim(strofreal(x/1e6, "%30." + strofreal(rd) + "f")) + "M")
    }
    if ((letter == "" & ax >= 1e3) | letter == "K") {
        return(strtrim(strofreal(x/1e3, "%30." + strofreal(rd) + "f")) + "K")
    }
    if (ax >= 1 | x == 0) {
        return(strtrim(strofreal(x, "%30." + strofreal(rd) + "f")))
    }
    fz = -floor(log10(ax))
    return(strtrim(strofreal(x, "%30." + strofreal(rd + fz - 1) + "f")))
}

string scalar _ci_letterof(string scalar s)
{
    string scalar c
    c = substr(s, strlen(s), 1)
    if (c == "B" | c == "M" | c == "K") {
        return(c)
    }
    return("")
}

string scalar _ci_prettypct(real scalar x)
{
    string scalar s
    s = strtrim(strofreal(100*x, "%30.0f"))
    if (100*x >= 0) {
        s = "+" + s
    }
    return(s + "%")
}

void _ci_report_strings(real scalar digits)
{
    real matrix   S
    string scalar L1, L2, a1, a2
    S  = st_matrix("e(summary)")
    a1 = _ci_prettynum(S[1,1], "", digits)
    a2 = _ci_prettynum(S[2,1], "", digits)
    L1 = _ci_letterof(a1)
    L2 = _ci_letterof(a2)
    st_local("r_act1",    a1)
    st_local("r_act2",    a2)
    st_local("r_pred1",   _ci_prettynum(S[1,2],  L1, 2))
    st_local("r_pred2",   _ci_prettynum(S[2,2],  L2, 2))
    st_local("r_predlo1", _ci_prettynum(S[1,3],  L1, digits))
    st_local("r_predlo2", _ci_prettynum(S[2,3],  L2, digits))
    st_local("r_predup1", _ci_prettynum(S[1,4],  L1, digits))
    st_local("r_predup2", _ci_prettynum(S[2,4],  L2, digits))
    st_local("r_abs1",    _ci_prettynum(S[1,6],  L1, digits))
    st_local("r_abs2",    _ci_prettynum(S[2,6],  L2, digits))
    st_local("r_abslo1",  _ci_prettynum(S[1,7],  L1, digits))
    st_local("r_absup1",  _ci_prettynum(S[1,8],  L1, digits))
    st_local("r_rel1",    _ci_prettypct(S[1,10]))
    st_local("r_rello1",  _ci_prettypct(S[1,11]))
    st_local("r_relup1",  _ci_prettypct(S[1,12]))
}

// -------------------------------------------------------------------------
// Small helpers
// -------------------------------------------------------------------------
real scalar _ci_firstge(real colvector tv, real scalar val)
{
    real scalar i
    for (i = 1; i <= rows(tv); i++) {
        if (tv[i] >= val) {
            return(i)
        }
    }
    return(.)
}

real scalar _ci_lastle(real colvector tv, real scalar val)
{
    real scalar i, k
    k = .
    for (i = 1; i <= rows(tv); i++) {
        if (tv[i] <= val) {
            k = i
        }
    }
    return(k)
}

real scalar _ci_illcond(real colvector y)
{
    real colvector nz
    nz = select(y, y :< .)
    if (length(nz) < 3) {
        return(1)
    }
    if (sqrt(variance(nz)) <= 0) {
        return(1)
    }
    return(0)
}

real scalar _ci_sdcol(real colvector x)
{
    real colvector nz
    real scalar    s
    nz = select(x, x :< .)
    if (length(nz) < 2) {
        return(1)
    }
    s = sqrt(variance(nz))
    if (s <= 0 | s >= .) {
        return(1)
    }
    return(s)
}

// -------------------------------------------------------------------------
//                              MAIN DRIVER
// -------------------------------------------------------------------------
void _ci_main()
{
    struct cissm scalar M

    external real matrix    _causalimpact_series
    external real colvector _causalimpact_rows

    real colvector idx, tv, yorig, ymod, keeprow, yfit, yobs, obsi, obsi2
    real colvector mu_c, sd_c, RQRd, yker, mucon, regcon, ydot, Xtyf
    real colvector bpr, pivec, logpi, log1mpi, gam, beta, sobs_d, slev_d
    real colvector ssea_d, sbet_d, pmean, plo, pup, cpred, cplo, cpup, cumy
    real colvector post_idx, ypost, rmean, rsum, dmean, dsum, relc
    real colvector pe, pel, peu, ce, cel, ceu, ae, rows_full, bmean
    real colvector Tadv1, a1v, advv, gj, bj, bnz, s2bet, dynsum, dynn
    real matrix    X, Xmod, Xd, Xobs, Zmat, alpha_t, statesam, ysam
    real matrix    betasam, gamsam, cumsam, XtXf, Om, SUM, BB, VV, INCL
    real matrix    bdev, Tadv, Thold, P1m
    string scalar  tousev, depv
    string rowvector covn, xnames

    real scalar n, i, j, t, m, nx, hasx, d, dss, Jd, ndraw, keep, it, kk
    real scalar pre1, pre2, post1, post2, i_pre1, i_pre2, i_post1, i_post2
    real scalar r_pre2, r_post1, r_post2, nobs, np
    real scalar niter, burnin, alpha, plevsd, S, dur, dynreg, stdz, addcon
    real scalar maxflips, msize, r2, pdf, ginfo, dshr, sigupper, dots, dotev
    real scalar ymu, ysd, sdy, s2obs, s2lev, s2seas, mj, sj, sxj
    real scalar dlev, sslev, uplev, dsea, sssea, upsea, dobs, ssobs, upobs
    real scalar dbet, upbet, nueps, seps, statreg, sumsq, nsum, eta, acc
    real scalar plow, pupp, run, lastc, cacc, col, acc2, sk, ysumobs
    real scalar ymeanobs, predsum_post, predmean_post, cge, cle, pval
    real scalar inpre, inpost, nrun, yty2, y1
    real colvector xj

    // ---------------------------------------------------------------- read
    tousev = st_local("touse")
    idx    = selectindex(st_data(., tousev) :!= 0)
    n      = length(idx)

    depv = st_local("depvar")
    covn = tokens(st_local("covars"))
    nx    = cols(covn)
    hasx = (nx > 0)

    tv    = st_data(idx, st_local("tvar"))
    yorig = st_data(idx, depv)
    X     = J(n, 0, 0)
    if (hasx) {
        X = st_data(idx, covn)
    }

    niter    = strtoreal(st_local("niter"))
    burnin   = strtoreal(st_local("burnin"))
    alpha    = strtoreal(st_local("alpha"))
    plevsd   = strtoreal(st_local("priorlevelsd"))
    S        = strtoreal(st_local("nseasons"))
    dur      = strtoreal(st_local("seasonduration"))
    dynreg   = strtoreal(st_local("dynreg"))
    stdz     = strtoreal(st_local("stdz"))
    addcon   = strtoreal(st_local("addcon"))
    maxflips = strtoreal(st_local("maxflips"))
    msize    = strtoreal(st_local("modelsize"))
    r2       = strtoreal(st_local("r2"))
    pdf      = strtoreal(st_local("priordf"))
    ginfo    = strtoreal(st_local("ginfo"))
    dshr     = strtoreal(st_local("dshrinkage"))
    sigupper = strtoreal(st_local("sigmaupper"))
    dots     = strtoreal(st_local("dots"))
    pre1     = strtoreal(st_local("pre1"))
    pre2     = strtoreal(st_local("pre2"))
    post1    = strtoreal(st_local("post1"))
    post2    = strtoreal(st_local("post2"))
    if (sigupper < 0) {
        sigupper = .
    }

    // ---------------------------------------------------------------- S1
    i_pre1  = _ci_firstge(tv, pre1)
    i_pre2  = _ci_lastle(tv, pre2)
    i_post1 = _ci_firstge(tv, post1)
    i_post2 = _ci_lastle(tv, post2)
    if (i_pre1 >= . | i_pre2 >= . | i_post1 >= . | i_post2 >= .) {
        st_local("e_ci_ok", "0")
        st_local("e_ci_msg", "each period must cover at least one data point")
        return
    }
    if (i_post2 > n) {
        i_post2 = n
    }

    for (i = 1; i <= n; i++) {
        if (yorig[i] < .) {
            break
        }
    }
    if (i > i_pre1) {
        i_pre1 = i
    }
    if (i_pre2 - i_pre1 < 2) {
        st_local("e_ci_ok", "0")
        st_local("e_ci_msg", "the pre-period must span at least 3 time points")
        return
    }
    if (i_post1 <= i_pre2) {
        st_local("e_ci_ok", "0")
        st_local("e_ci_msg", "the post-period must start after the pre-period ends")
        return
    }

    // ---------------------------------------------------------------- S2
    keeprow = (i_pre1::n)
    m       = length(keeprow)
    ymod    = yorig[keeprow]
    Xmod    = J(m, 0, 0)
    if (hasx) {
        Xmod = X[keeprow, .]
    }
    r_pre2  = i_pre2  - i_pre1 + 1
    r_post1 = i_post1 - i_pre1 + 1
    r_post2 = i_post2 - i_pre1 + 1

    if (_ci_illcond(ymod[1::r_pre2])) {
        st_local("e_ci_ok", "0")
        st_local("e_ci_msg", "the response is constant, all missing, or has fewer than 3 non-missing values in the pre-period")
        return
    }

    ymu  = 0
    ysd  = 1
    mu_c = J(max((nx,1)), 1, 0)
    sd_c = J(max((nx,1)), 1, 1)
    if (stdz) {
        ymu = mean(select(ymod[1::r_pre2], ymod[1::r_pre2] :< .))
        ysd = _ci_sdcol(ymod[1::r_pre2])
        ymod = (ymod :- ymu)/ysd
        if (hasx) {
            for (j = 1; j <= nx; j++) {
                xj = Xmod[., j]
                mj = mean(select(xj[1::r_pre2], xj[1::r_pre2] :< .))
                sj = _ci_sdcol(xj[1::r_pre2])
                Xmod[., j] = (xj :- mj)/sj
                mu_c[j] = mj
                sd_c[j] = sj
            }
        }
    }

    // Design matrix; bsts includes an intercept in the regression component.
    Xd     = Xmod
    xnames = covn
    if (hasx & addcon) {
        Xd     = (Xmod, J(m, 1, 1))
        xnames = (covn, "_cons")
    }
    Jd = cols(Xd)

    // ---------------------------------------------------------------- S3
    yfit = ymod
    for (t = r_pre2 + 1; t <= m; t++) {
        yfit[t] = .
    }
    nobs = 0
    for (t = 1; t <= m; t++) {
        if (yfit[t] < .) {
            nobs = nobs + 1
        }
    }

    sdy = _ci_sdcol(yfit)

    // ---------------------------------------------------------------- S4-S8
    dss = 1
    if (S > 1) {
        dss = S
    }
    d = dss
    if (dynreg) {
        d = dss + Jd
    }

    Tadv  = J(d, d, 0)
    Thold = J(d, d, 0)
    Tadv[1,1]  = 1
    Thold[1,1] = 1
    if (S > 1) {
        for (i = 2; i <= dss; i++) {
            Tadv[2,i] = -1
        }
        for (i = 3; i <= dss; i++) {
            Tadv[i,i-1] = 1
        }
        for (i = 2; i <= dss; i++) {
            Thold[i,i] = 1
        }
    }
    if (dynreg) {
        for (i = dss + 1; i <= d; i++) {
            Tadv[i,i]  = 1
            Thold[i,i] = 1
        }
    }

    advv = J(m, 1, 1)
    if (dur > 1) {
        for (t = 1; t <= m; t++) {
            advv[t] = (mod(t, dur) == 0)
        }
    }

    Zmat = J(d, m, 0)
    for (t = 1; t <= m; t++) {
        Zmat[1,t] = 1
        if (S > 1) {
            Zmat[2,t] = 1
        }
        if (dynreg) {
            for (j = 1; j <= Jd; j++) {
                Zmat[dss + j, t] = Xd[t, j]
            }
        }
    }

    y1 = 0
    for (t = 1; t <= m; t++) {
        if (yfit[t] < .) {
            y1 = yfit[t]
            break
        }
    }
    a1v      = J(d, 1, 0)
    a1v[1]   = y1
    P1m      = I(d)*(sdy*sdy)
    s2bet    = J(max((Jd,1)), 1, 0)
    if (dynreg) {
        for (j = 1; j <= Jd; j++) {
            sxj = _ci_sdcol(Xd[., j])
            P1m[dss + j, dss + j] = (sdy/sxj)^2
            s2bet[j] = (0.01*sdy/sxj)^2
        }
    }

    M.Tadv  = Tadv
    M.Thold = Thold
    M.a1    = a1v
    M.P1    = P1m
    M.adv   = advv
    M.d     = d
    M.S     = S
    M.dur   = dur
    M.m     = m
    M.nx     = Jd
    M.dyn   = dynreg
    M.dss   = dss

    // ---------------------------------------------------------------- S9
    dlev  = 32
    sslev = dlev*(plevsd*sdy)^2
    uplev = sdy

    dsea  = 0.01
    sssea = dsea*(0.01*sdy)^2
    upsea = sdy

    dbet  = 32
    upbet = sdy

    if (dynreg) {
        dobs  = 32
        ssobs = dobs*(plevsd*sdy)^2
        upobs = 0.1*sdy
    }
    else {
        dobs  = 0.01
        ssobs = dobs*(sdy*sdy)
        upobs = 1.2*sdy
    }

    // ---------------------------------------------------------------- S10-S11
    statreg = (hasx & !dynreg)
    gam     = J(max((Jd,1)), 1, 0)
    beta    = J(max((Jd,1)), 1, 0)
    XtXf    = J(1, 1, 0)
    Om      = J(1, 1, 0)
    bpr     = J(max((Jd,1)), 1, 0)
    pivec   = J(max((Jd,1)), 1, 0)
    logpi   = J(max((Jd,1)), 1, 0)
    log1mpi = J(max((Jd,1)), 1, 0)
    nueps   = pdf
    seps    = pdf*(1 - r2)*sdy*sdy

    if (statreg) {
        obsi = selectindex(yfit :< .)
        Xobs = Xd[obsi, .]
        XtXf = quadcross(Xobs, Xobs)
        Om   = (ginfo/nobs)*((1 - dshr)*XtXf + dshr*diag(diagonal(XtXf)))
        for (j = 1; j <= Jd; j++) {
            if (Om[j,j] <= 0) {
                Om[j,j] = 1e-8
            }
        }
        bpr     = J(Jd, 1, 0)
        pivec   = J(Jd, 1, min((1, msize/Jd)))
        logpi   = log(pivec)
        log1mpi = log(1 :- pivec)
        gam     = J(Jd, 1, 1)
        beta    = J(Jd, 1, 0)
    }

    // ---------------------------------------------------------------- init
    s2obs  = (0.1*sdy)^2
    s2lev  = (plevsd*sdy)^2
    s2seas = (0.01*sdy)^2

    ndraw = niter - burnin
    if (ndraw*m > 60000000) {
        st_local("e_ci_ok", "0")
        st_local("e_ci_msg", "the requested MCMC sample is too large for memory; reduce niter() or the sample")
        return
    }

    statesam = J(ndraw, m, 0)
    ysam     = J(ndraw, m, 0)
    betasam  = J(ndraw, max((Jd,1)), 0)
    gamsam   = J(ndraw, max((Jd,1)), 0)
    sobs_d   = J(ndraw, 1, 0)
    slev_d   = J(ndraw, 1, 0)
    ssea_d   = J(ndraw, 1, 0)

    dotev = ceil(niter/50)
    if (dots) {
        printf("{txt}Gibbs sampler: ")
        displayflush()
    }

    keep = 0
    for (it = 1; it <= niter; it++) {

        // ------------------------------------------------------------ S12
        RQRd    = J(d, 1, 0)
        RQRd[1] = s2lev
        if (S > 1) {
            RQRd[2] = s2seas
        }
        if (dynreg) {
            for (j = 1; j <= Jd; j++) {
                RQRd[dss + j] = s2bet[j]
            }
        }

        yker = yfit
        if (statreg) {
            for (t = 1; t <= m; t++) {
                if (yfit[t] < .) {
                    yker[t] = yfit[t] - Xd[t, .]*beta
                }
            }
        }

        alpha_t = _ci_simsmooth(yker, M, Zmat, s2obs, RQRd)

        // ------------------------------------------------------------ S9
        sumsq = 0
        nsum  = 0
        for (t = 1; t <= m - 1; t++) {
            eta   = alpha_t[1, t+1] - alpha_t[1, t]
            sumsq = sumsq + eta*eta
            nsum  = nsum + 1
        }
        s2lev = _ci_draw_sig2(dlev, sslev, nsum, sumsq, uplev)
        if (s2lev >= . | s2lev <= 0) {
            s2lev = (plevsd*sdy)^2
        }

        if (S > 1) {
            sumsq = 0
            nsum  = 0
            for (t = 1; t <= m - 1; t++) {
                if (advv[t] == 1) {
                    eta = alpha_t[2, t+1]
                    for (i = 2; i <= dss; i++) {
                        eta = eta + alpha_t[i, t]
                    }
                    sumsq = sumsq + eta*eta
                    nsum  = nsum + 1
                }
            }
            s2seas = _ci_draw_sig2(dsea, sssea, nsum, sumsq, upsea)
            if (s2seas >= . | s2seas <= 0) {
                s2seas = (0.01*sdy)^2
            }
        }

        if (dynreg) {
            for (j = 1; j <= Jd; j++) {
                sumsq = 0
                nsum  = 0
                for (t = 1; t <= m - 1; t++) {
                    eta   = alpha_t[dss + j, t+1] - alpha_t[dss + j, t]
                    sumsq = sumsq + eta*eta
                    nsum  = nsum + 1
                }
                sxj = _ci_sdcol(Xd[., j])
                s2bet[j] = _ci_draw_sig2(dbet, dbet*(0.01*sdy/sxj)^2, nsum,
                                         sumsq, upbet/sxj)
                if (s2bet[j] >= . | s2bet[j] <= 0) {
                    s2bet[j] = (0.01*sdy/sxj)^2
                }
            }
        }

        // ------------------------------------------------------------ S13
        mucon = J(m, 1, 0)
        for (t = 1; t <= m; t++) {
            mucon[t] = alpha_t[1, t]
            if (S > 1) {
                mucon[t] = mucon[t] + alpha_t[2, t]
            }
        }
        regcon = J(m, 1, 0)
        if (dynreg) {
            for (t = 1; t <= m; t++) {
                acc = 0
                for (j = 1; j <= Jd; j++) {
                    acc = acc + Xd[t, j]*alpha_t[dss + j, t]
                }
                regcon[t] = acc
            }
        }

        if (statreg) {
            ydot = J(m, 1, .)
            for (t = 1; t <= m; t++) {
                if (yfit[t] < .) {
                    ydot[t] = yfit[t] - mucon[t]
                }
            }
            obsi2 = selectindex(ydot :< .)
            Xtyf  = quadcross(Xd[obsi2, .], ydot[obsi2])
            yty2  = quadcross(ydot[obsi2], ydot[obsi2])
            _ci_ssvs(gam, beta, s2obs, XtXf, Xtyf, yty2, Om, bpr, nueps, seps,
                     nobs, logpi, log1mpi, pivec, maxflips, sigupper)
            regcon = Xd*beta
        }
        else {
            sumsq = 0
            nsum  = 0
            for (t = 1; t <= m; t++) {
                if (yfit[t] < .) {
                    eta   = yfit[t] - mucon[t] - regcon[t]
                    sumsq = sumsq + eta*eta
                    nsum  = nsum + 1
                }
            }
            s2obs = _ci_draw_sig2(dobs, ssobs, nsum, sumsq, upobs)
            if (s2obs >= . | s2obs <= 0) {
                s2obs = sumsq/max((nsum, 1))
            }
        }

        // ------------------------------------------------------------ S14
        if (it > burnin) {
            keep = keep + 1
            for (t = 1; t <= m; t++) {
                statesam[keep, t] = mucon[t] + regcon[t]
            }
            ysam[keep, .] = statesam[keep, .] + rnormal(1, m, 0, sqrt(s2obs))
            if (Jd > 0) {
                if (statreg) {
                    betasam[keep, .] = beta'
                    gamsam[keep, .]  = gam'
                }
                if (dynreg) {
                    for (j = 1; j <= Jd; j++) {
                        betasam[keep, j] = mean(alpha_t[dss + j, .]')
                        gamsam[keep, j]  = 1
                    }
                }
            }
            sobs_d[keep] = sqrt(s2obs)
            slev_d[keep] = sqrt(s2lev)
            ssea_d[keep] = sqrt(s2seas)
        }

        if (dots) {
            if (mod(it, dotev) == 0) {
                printf(".")
                displayflush()
            }
        }
    }
    if (dots) {
        printf("{txt} done\n")
        displayflush()
    }

    // -------------------------------------------------- undo S2 before S16
    if (stdz) {
        statesam = statesam*ysd :+ ymu
        ysam     = ysam*ysd :+ ymu
    }
    yobs = yorig[keeprow]

    // ---------------------------------------------------------------- S14
    plow  = alpha/2
    pupp  = 1 - alpha/2
    pmean = J(m, 1, 0)
    plo   = J(m, 1, 0)
    pup   = J(m, 1, 0)
    for (t = 1; t <= m; t++) {
        pmean[t] = mean(statesam[., t])
        plo[t]   = _ci_qcol(ysam[., t], plow)
        pup[t]   = _ci_qcol(ysam[., t], pupp)
    }

    // ---------------------------------------------------------------- S16
    cpred = J(m, 1, .)
    cplo  = J(m, 1, .)
    cpup  = J(m, 1, .)
    cumy  = J(m, 1, .)
    run   = 0
    lastc = 0
    for (t = 1; t <= m; t++) {
        if (yobs[t] < .) {
            run     = run + yobs[t]
            cumy[t] = run
        }
        if (t < r_post1) {
            if (yobs[t] < .) {
                cpred[t] = run
                cplo[t]  = run
                cpup[t]  = run
                lastc    = run
            }
        }
    }
    cumsam = J(ndraw, m - r_post1 + 1, 0)
    for (kk = 1; kk <= ndraw; kk++) {
        acc2 = 0
        col  = 0
        for (t = r_post1; t <= m; t++) {
            col             = col + 1
            acc2            = acc2 + ysam[kk, t]
            cumsam[kk, col] = acc2 + lastc
        }
    }
    cacc = lastc
    col  = 0
    for (t = r_post1; t <= m; t++) {
        col      = col + 1
        cacc     = cacc + pmean[t]
        cpred[t] = cacc
        cplo[t]  = _ci_qcol(cumsam[., col], plow)
        cpup[t]  = _ci_qcol(cumsam[., col], pupp)
    }

    // ---------------------------------------------------------------- S18
    post_idx = (r_post1::r_post2)
    ypost    = yobs[post_idx]
    np       = length(post_idx)
    ysumobs  = quadsum(ypost)
    ymeanobs = ysumobs/np

    rmean = J(ndraw, 1, 0)
    rsum  = J(ndraw, 1, 0)
    dmean = J(ndraw, 1, 0)
    dsum  = J(ndraw, 1, 0)
    relc  = J(ndraw, 1, 0)
    for (kk = 1; kk <= ndraw; kk++) {
        sk        = quadsum(ysam[kk, post_idx])
        rsum[kk]  = sk
        rmean[kk] = sk/np
        dsum[kk]  = ysumobs - sk
        dmean[kk] = ymeanobs - sk/np
        relc[kk]  = ysumobs/sk - 1
    }
    predsum_post  = quadsum(pmean[post_idx])
    predmean_post = predsum_post/np

    SUM = J(2, 15, 0)
    SUM[1,1]  = ymeanobs
    SUM[2,1]  = ysumobs
    SUM[1,2]  = predmean_post
    SUM[2,2]  = predsum_post
    SUM[1,3]  = _ci_qcol(rmean, plow)
    SUM[2,3]  = _ci_qcol(rsum,  plow)
    SUM[1,4]  = _ci_qcol(rmean, pupp)
    SUM[2,4]  = _ci_qcol(rsum,  pupp)
    SUM[1,5]  = sqrt(variance(rmean))
    SUM[2,5]  = sqrt(variance(rsum))
    SUM[1,6]  = ymeanobs - predmean_post
    SUM[2,6]  = ysumobs  - predsum_post
    SUM[1,7]  = _ci_qcol(dmean, plow)
    SUM[2,7]  = _ci_qcol(dsum,  plow)
    SUM[1,8]  = _ci_qcol(dmean, pupp)
    SUM[2,8]  = _ci_qcol(dsum,  pupp)
    SUM[1,9]  = sqrt(variance(dmean))
    SUM[2,9]  = sqrt(variance(dsum))
    // R forces the average and cumulative relative effect to be identical.
    SUM[1,10] = mean(relc)
    SUM[2,10] = mean(relc)
    SUM[1,11] = _ci_qcol(relc, plow)
    SUM[2,11] = _ci_qcol(relc, plow)
    SUM[1,12] = _ci_qcol(relc, pupp)
    SUM[2,12] = _ci_qcol(relc, pupp)
    SUM[1,13] = sqrt(variance(relc))
    SUM[2,13] = sqrt(variance(relc))
    SUM[1,14] = alpha
    SUM[2,14] = alpha

    cge = 1
    cle = 1
    for (kk = 1; kk <= ndraw; kk++) {
        if (rsum[kk] >= ysumobs) {
            cge = cge + 1
        }
        if (rsum[kk] <= ysumobs) {
            cle = cle + 1
        }
    }
    pval      = min((cge, cle))/(ndraw + 1)
    SUM[1,15] = pval
    SUM[2,15] = pval

    st_matrix(st_local("SUM"), SUM)
    st_matrixrowstripe(st_local("SUM"), (J(2,1,""), ("Average" \ "Cumulative")))
    st_matrixcolstripe(st_local("SUM"), (J(15,1,""),
        ("Actual" \ "Pred" \ "Pred_lower" \ "Pred_upper" \ "Pred_sd" \
         "AbsEffect" \ "AbsEffect_lower" \ "AbsEffect_upper" \ "AbsEffect_sd" \
         "RelEffect" \ "RelEffect_lower" \ "RelEffect_upper" \ "RelEffect_sd" \
         "alpha" \ "p")))

    // ------------------------------------------------------------ e(b), e(V)
    if (Jd > 0) {
        bmean = J(Jd, 1, 0)
        for (j = 1; j <= Jd; j++) {
            bmean[j] = mean(betasam[., j])
        }
        bdev = betasam :- bmean'
        BB   = bmean'
        VV   = quadcross(bdev, bdev)/(ndraw - 1)

        INCL = J(Jd, 4, 0)
        for (j = 1; j <= Jd; j++) {
            gj = gamsam[., j]
            bj = betasam[., j]
            INCL[j,1] = mean(gj)
            INCL[j,2] = mean(bj)
            INCL[j,3] = sqrt(variance(bj))
            bnz = select(bj, gj :> 0)
            if (length(bnz) > 0) {
                INCL[j,4] = mean(bnz :> 0)
            }
            else {
                INCL[j,4] = .
            }
        }
        st_matrix(st_local("BB"), BB)
        st_matrix(st_local("VV"), VV)
        st_matrixcolstripe(st_local("BB"), (J(Jd,1,""), xnames'))
        st_matrixcolstripe(st_local("VV"), (J(Jd,1,""), xnames'))
        st_matrixrowstripe(st_local("VV"), (J(Jd,1,""), xnames'))
        st_matrix(st_local("INCL"), INCL)
        st_matrixrowstripe(st_local("INCL"), (J(Jd,1,""), xnames'))
        st_matrixcolstripe(st_local("INCL"), (J(4,1,""),
            ("P_include" \ "PostMean" \ "PostSD" \ "P_positive")))
    }

    // ---------------------------------------------------------------- S15-S17
    rows_full = idx[keeprow]

    st_store(rows_full, st_local("s_resp"),    yobs)
    st_store(rows_full, st_local("s_cumresp"), cumy)
    st_store(rows_full, st_local("s_pred"),    pmean)
    st_store(rows_full, st_local("s_predlo"),  plo)
    st_store(rows_full, st_local("s_predup"),  pup)
    st_store(rows_full, st_local("s_cpred"),   cpred)
    st_store(rows_full, st_local("s_cpredlo"), cplo)
    st_store(rows_full, st_local("s_cpredup"), cpup)

    pe   = J(m, 1, .)
    pel  = J(m, 1, .)
    peu  = J(m, 1, .)
    ce   = J(m, 1, .)
    cel  = J(m, 1, .)
    ceu  = J(m, 1, .)
    ae   = J(m, 1, .)
    nrun = 0
    for (t = 1; t <= m; t++) {
        inpre  = (t <= r_pre2)
        inpost = (t >= r_post1 & t <= r_post2)
        if (inpre | inpost) {
            if (yobs[t] < .) {
                pe[t]  = yobs[t] - pmean[t]              // S15  eq. (2.15)
                pel[t] = yobs[t] - pup[t]
                peu[t] = yobs[t] - plo[t]
            }
            if (cumy[t] < . & cpred[t] < .) {            // S16  eq. (2.16)
                ce[t]  = cumy[t] - cpred[t]
                cel[t] = cumy[t] - cpup[t]
                ceu[t] = cumy[t] - cplo[t]
            }
        }
        if (inpost) {                                    // S17  eq. (2.17)
            nrun = nrun + 1
            if (ce[t] < .) {
                ae[t] = ce[t]/nrun
            }
        }
    }

    st_store(rows_full, st_local("s_eff"),    pe)
    st_store(rows_full, st_local("s_efflo"),  pel)
    st_store(rows_full, st_local("s_effup"),  peu)
    st_store(rows_full, st_local("s_ceff"),   ce)
    st_store(rows_full, st_local("s_cefflo"), cel)
    st_store(rows_full, st_local("s_ceffup"), ceu)
    st_store(rows_full, st_local("s_avgeff"), ae)

    // Keep the series available to causalimpact_p (predict).
    _causalimpact_series = (pmean, plo, pup, cpred, cplo, cpup, pe, pel, peu,
                            ce, cel, ceu, yobs, cumy, ae)
    _causalimpact_rows   = rows_full

    // ---------------------------------------------------------------- out
    st_local("e_ci_N",      strofreal(n))
    st_local("e_ci_npre",   strofreal(i_pre2 - i_pre1 + 1))
    st_local("e_ci_npost",  strofreal(np))
    st_local("e_ci_nmodel", strofreal(m))
    st_local("e_ci_p",      strofreal(pval,      "%20.15g"))
    st_local("e_ci_act1",   strofreal(SUM[1,1],  "%20.15g"))
    st_local("e_ci_act2",   strofreal(SUM[2,1],  "%20.15g"))
    st_local("e_ci_pred1",  strofreal(SUM[1,2],  "%20.15g"))
    st_local("e_ci_pred2",  strofreal(SUM[2,2],  "%20.15g"))
    st_local("e_ci_abs1",   strofreal(SUM[1,6],  "%20.15g"))
    st_local("e_ci_abs2",   strofreal(SUM[2,6],  "%20.15g"))
    st_local("e_ci_abssd1", strofreal(SUM[1,9],  "%20.15g"))
    st_local("e_ci_abssd2", strofreal(SUM[2,9],  "%20.15g"))
    st_local("e_ci_rel1",   strofreal(SUM[1,10], "%20.15g"))
    st_local("e_ci_relsd1", strofreal(SUM[1,13], "%20.15g"))
    st_local("e_ci_sobs",   strofreal(mean(sobs_d)*ysd, "%20.15g"))
    st_local("e_ci_slev",   strofreal(mean(slev_d)*ysd, "%20.15g"))
    st_local("e_ci_sseas",  strofreal(mean(ssea_d)*ysd, "%20.15g"))
}

end