hospitalization.do

// Requires egenmore package, which can be installed with:
// ssc install egenmore

import delimited https://covidtracking.com/api/v1/us/daily.csv, clear

gen newdate = date(strofreal(date,"%8.0f"),"YMD")
format newdate %tdCCYY-NN-DD
drop date
rename newdate date
label var date "Date"

tsset date
label var hospitalizedcurrently "Currently hospitalized"

gen dailypositives = d.positive
gen positiverate = positive/(positive+negative) // exclude pending tests in calculating positiverate
gen estnewinfections = max(dailypositives*20*positiverate,dailypositives) // assume positive rate that implies full test coverage = 0.05

tsappend, add(45) // add 45 days to time series

set more off // this block does some forecasting
quietly sum date if estnewinfections != .
local latest = r(max)
quietly reg estnewinfections date if date > `latest' - 17 // use last 2.5 weeks of data to predict infections going forward
estimates store recentinfections
forecast create simplemodel, replace
forecast estimates recentinfections
forecast solve, begin(`latest') end(`latest'+45)

tssmooth ma estmainfections = f_estnewinfections, window(14 1) // create moving average of last 14 days + today
gen estinfectionslast15days = 15*estmainfections // multiply by 15 to generate the number of infections in the last 15 days

gen estcuminfections = sum(estnewinfections) // this block estimates the hospitalization rate for the population as a whole
gen pcthospitalized = 100*hospitalizedcumulative/estcuminfections
egen esthospitalizationrate = lastnm(pcthospitalized) // uses the latest cumulative hospitalization rate
egen temp = mean(esthospitalizationrate)
replace esthospitalizationrate = temp/100
drop temp // needed to use this temp variable to extend hospitalization rate to appended days

gen predictedinhospital = esthospitalizationrate*l1.estinfectionslast15days
label var predictedinhospital "Predicted in hospital"

tsline hospitalizedcurrently predictedinhospital, xsize(16) ysize(9) scale(0.6) ylabel(,format(%15.0fc)) title("Current and predicted COVID-19 hospitalizations in the US") xtitle("") note("A very simple model by @elidourado. Extrapolates infections by last 2.5 weeks of data. Raw data from covidtracking.com." "Stata code available at https://gist.github.com/elidourado/9b89c71e9051a97a5d53209c357d200d.")
graph export hospitalization-prediction.png, as(png) width(960) height(540) replace