* Script: Cornoavirus.do * Author: Anna Ziff (anna.ziff@duke.edu) * Original date: February 10, 2020 * Purpose: This script analyzes the WHO reports on neo-CoV clear all * Open log file cd $results cap log close log using "CoronavirusLog.log", replace text * Set seed set seed 1024 * Environment global projects : env projects global repo = "${projects}/CoV" global scripts = "${repo}/Scripts" global results = "${repo}/Results" global data = "${repo}/Data" * Import data cd $data import excel using "coronavirus", cellrange(G3:P30) firstrow cd $results * Clean data rename K sumconf rename M chinadeaths rename O totaldeaths rename P newdeaths gen report2 = report^2 gen report3 = report^3 gen chinarate = chinadeaths/chinaconf gen totalrate = totaldeaths/totalconf gen logtotalconf = log(totalconf) gen logreport = log(report) gen logtotaldeaths = log(totaldeaths) lab var date "Date" lab var report "Time" lab var chinaconf "Confirmed Cases: China" lab var nonchinaconf "Confirmed Cases: Non-China" lab var sumconf "Summed Confirmed Cases (calculated)" lab var totalconf "Total Confirmed Cases" lab var chinadeaths "Deaths: China" lab var nonchinadeaths "Deaths: Non-China" lab var totaldeaths "Total Deaths" lab var newdeaths "New Deaths (calculated)" lab var report2 "Time^2" lab var chinarate "Death Rate: China" lab var totalrate "Total Death Rate" lab var logtotalconf "log(Total Confirmed Cases)" lab var logreport "log(t)" lab var logtotaldeaths "log(Total Deaths)" format date %tdnn/dd/YY ************** * Regression * ************** * Quadratic Fit reg totaldeaths report report2 if report != 24, robust eststo quad esttab quad using DeathsQuadratic.tex, replace booktabs label nomtitles b se r2 * Test power-law fit omitting more t forvalues t = 1/27 { if `t' <= 21 { reg logtotaldeaths logreport if report >= `t' & report != 24, robust // Save results matrix tmpb = e(b) matrix tmpv = e(V) local b`t' = tmpb[1,1] local v`t' = sqrt(tmpv[1,1]) if `t' == 1 { matrix b = `b`t'' matrix v = `v`t'' } else { matrix b = (b \ `b`t'') matrix v = (v \ `v`t'') } } if `t' >= 11 { reg logtotaldeaths logreport if report <= `t' & report >= 8 & report != 24, robust // Save results matrix tmpb = e(b) matrix tmpv = e(V) local b`t' = tmpb[1,1] local v`t' = sqrt(tmpv[1,1]) if `t' == 11 { matrix b_back = `b`t'' matrix v_back = `v`t'' } else { matrix b_back = (b_back \ `b`t'') matrix v_back = (v_back \ `v`t'') } } } matrix powercoef = (b , v) // Matrix of power-law coefficients and s.e. matrix powercoef_back = (b_back, v_back) * Nonparametric: Test power law coefficient omitting the first time period // Get point estimate qui reg logtotaldeaths logreport if report >= 2 & report != 24 matrix tmpb = e(b) matrix b2 = tmpb[1,1] qui reg logtotaldeaths logreport if report >= 8 & report != 24 matrix tmpb = e(b) matrix b6 = tmpb[1,1] qui reg totaldeaths report report2 if report != 24 matrix tmpb = e(b) matrix quad1 = tmpb[1,1] matrix quad2 = tmpb[1,2] matrix quad3 = tmpb[1,3] // Bootstrap (set bootstrap number here) forvalues b = 1/5000 { preserve bsample // Randomly sample with replacement // Power-law: Regress for periods 2-25 qui reg logtotaldeaths logreport if report >= 2 & report != 24 matrix tmpb = e(b) matrix b2 = (b2 \ tmpb[1,1]) // Power-law: Regress for periods 8-25 qui reg logtotaldeaths logreport if report >= 8 & report != 24 matrix tmpb = e(b) matrix b6 = (b6 \ tmpb[1,1]) // Quadratic qui reg totaldeaths report report2 if report != 24 matrix quadb = e(b) matrix quad1 = (quad1 \ quadb[1,1]) matrix quad2 = (quad2 \ quadb[1,2]) matrix quad3 = (quad3 \ quadb[1,3]) // Print update every 100 bootstraps - useful for long runs if mod(`b',100) == 0 di "Bootstrap: `b'" restore } matrix all = (b2, b6, quad1, quad2, quad3) preserve clear svmat all gen n = _n // Summarize bootstraps omitting the point estimate (n = 1) // Power-law exponent for periods 2-26 sum all1 if n > 1, detail // Power-law exponent for periods 8-26 sum all2 if n > 1, detail // Quadratic coefficient for exopnential fit sum all3 if n > 1, detail // Linear coefficient for exponential fit sum all4 if n > 1, detail // Constant for exponential fit sum all5 if n > 1, detail restore ********* * Graph * ********* * Raw trend of confirmed cases # delimit ; twoway (connected totalconf date, lcol(black) mcol(black)), graphregion(color(white)) ylabel(, glcolor(gs13) angle(0)); # delimit cr graph export "ConfirmedTotal.eps", replace * Raw trend of deaths # delimit ; twoway (connected totaldeaths date, lcol(black) mcol(black)), graphregion(color(white)) ylabel(, glcolor(gs13) angle(0)); # delimit cr graph export "DeathsTotal.eps", replace * Deaths vs. time on log-log qui reg logtotaldeaths logreport if report >= 8 & report != 24 predict raw gen exp = exp(raw) # delimit ; twoway (line exp report if report >= 8, lcol(red) lwidth(1.1)) (connected totaldeaths report if report != 24, lcol(black) mcol(black) msymb(Oh)) , yscale(log) xscale(log) graphregion(color(white)) ylabel(, glcolor(gs13) angle(0)) xlabel(, grid glcolor(gs13)) legend(order(2 1) label(1 "Power-law Fit") label(2 "Total Deaths")); # delimit cr graph export "DeathsTotal_loglog.eps", replace * Log deaths vs. time # delimit ; twoway (connected logtotaldeaths date if report != 24, lcol(black) mcol(black)), graphregion(color(white)) ylabel(, glcolor(gs13) angle(0)); # delimit cr graph export "DeathsTotal_log.eps", replace * Graph the power coefficient over N clear svmat powercoef gen t = _n gen cilower = powercoef1 - invttail(26-t+1-2,0.05)*powercoef2 gen ciupper = powercoef1 + invttail(26-t+1-2,0.05)*powercoef2 # delimit ; twoway (rcap cilower ciupper t, lcol(black)) (scatter powercoef1 t, mcol(black) lcol(black)), graphregion(color(white)) ytitle("Power-law Exponent") xtitle("Periods Analyzed") xlabel(1 "1-27" 2 "2-27" 3 "3-27" 4 "4-27" 5 "5-27" 6 "6-27" 7 "7-27" 8 "8-27" 9 "9-27" 10 "10-27" 11 "11-27" 12 "12-27" 13 "13-27" 14 "14-27" 15 "15-27" 16 "16-27" 17 "17-27" 18 "18-27" 19 "19-27" 20 "20-27" 21 "21-27", angle(30) labsize(*.8) grid glcolor(gs13)) ylabel(, glcolor(gs13) angle(0)) legend(label(1 "95% Confidence Interval") label(2 "Exponent")); # delimit cr graph export "PowerLawCoefficient.eps", replace clear svmat powercoef_back gen t = _n gen cilower = powercoef_back1 - invttail(t,0.05)*powercoef_back2 gen ciupper = powercoef_back1 + invttail(t,0.05)*powercoef_back2 # delimit ; twoway (rcap cilower ciupper t if t != 14, lcol(black)) (scatter powercoef_back1 t if t != 14, mcol(black) lcol(black)), graphregion(color(white)) ytitle("Power-law Exponent") xtitle("Periods Analyzed") xlabel(1 "8-11" 2 "8-12" 3 "8-13" 4 "8-14" 5 "8-15" 6 "8-16" 7 "8-17" 8 "8-18" 9 "8-19" 10 "8-20" 11 "8-21" 12 "8-22" 13 "8-23" 15 "8-25" 16 "8-26" 17 "8-27", angle(30) labsize(*.8) grid glcolor(gs13)) ylabel(, glcolor(gs13) angle(0)) legend(label(1 "95% Confidence Interval") label(2 "Exponent")); # delimit cr graph export "PowerLawCoefficient_back.eps", replace * Close log log close