'Simulation of SEIR model via add-in seirmodel(pop="1500", sim="250", r0="2.5", sigma="0.15", gamma="0.10", day0="12/10/2019", susceptible="1490", exposed="10", infected="0", recovered="0", chart, results) 'Fitting of nonlinear curves to simulated data of infected people series cum_infected = @cumsum(infected) 'cumulated series c = 0.1 'set starting values equation eq_glogistic.ls(showopts) cum_infected = c(1)/(1+(c(2)/@abs(c(3)))*@exp(-c(4)*@trend))^@abs(c(3)) 'for generalized logistic c= 0.1 'reset the starting values equation eq_gompertz.ls(showopts) cum_infected = c(1)*@exp(-c(2)*@exp(-c(4)*@trend)) 'for Gompertz 'Estimation of observational models infected.tsepigrowth(accum="1", growth="1", model="1") 'for generalized logistic infected.tsepigrowth(accum="1", growth="1", model="2") 'for Gompertz 'Get Covid19 data from WHO database wfopen https://covid19.who.int/WHO-COVID-19-global-data.csv @freq D7 @id country @date(date_reported) smpl @all if country_code="TR" 'restrict the sample to Turkey %first = @otod(@ifirst(cumulative_cases)) %last = "7/12/2020" '@otod(@ilast(cumulative_cases)) 'end date is fixed to replicate results in the blog post %end = "12/31/2020" pagecreate(page=Turkey) d7 {%first} {%last} 'create a time series page copy(smpl="@all if country_code=""TR""", c=none) Untitled\cumulative_cases total_cases @src @date @dest @date 'copy cumulative cases pagestruct(end={%end}) 'extend the range to cover for forecasting period 'Fit Gompertz curve equation gompertz_eq.ls(showopts, numericderiv) total_cases = c(1)*@exp(-(@abs(c(2)))*@exp(-c(4)*@trend)) gompertz_eq.fit(e, g) Gompertz_curve_f 'make forecasts 'Estimate observational counterparts of Gompertz tsepigrowth(ser="total_cases", growth="1", model="2", sample="3/12/2020 12/31/2020") 'for Gompertz tsepigrowth(ser="total_cases", growth="1", model="2", sample="3/12/2020 12/31/2020", policy="6/1/2020 7/12/2020") 'for Gompertz with slope intervention tsepigrowth(ser="total_cases", growth="1", model="3", sample="3/12/2020 12/31/2020") 'for Dynamic Gompertz 'Save the smoothed signal series ssmodel.makesignals(t=smooth) Gompertz_obs_g ssmodel01.makesignals(t=smooth) Gompertz_obs_intervention_g ssmodel02.makesignals(t=smooth) Gompertz_obs_dynamic_g 'Visual comparison of fitted values of growth series growth = log(d(total_cases)/total_cases(-1)) series Gompertz_curve_g = log(d(Gompertz_curve_f)/Gompertz_curve_f(-1)) smpl %first %last graph grfit.line growth Gompertz_curve_g Gompertz_obs_g Gompertz_obs_intervention_g Gompertz_obs_dynamic_g 'Generate forecast series smpl %last %end series Gompertz_obs_f = @recode(total_cases<>NA,total_cases(-1)*(1+@exp(Gompertz_obs_g)),Gompertz_obs_f(-1)*(1+@exp(Gompertz_obs_g))) series Gompertz_obs_intervention_f = @recode(total_cases<>NA,total_cases(-1)*(1+@exp(Gompertz_obs_intervention_g)),Gompertz_obs_intervention_f(-1)*(1+@exp(Gompertz_obs_intervention_g))) series Gompertz_obs_dynamic_f = @recode(total_cases<>NA,total_cases(-1)*(1+@exp(Gompertz_obs_dynamic_g)),Gompertz_obs_dynamic_f(-1)*(1+@exp(Gompertz_obs_dynamic_g))) 'Visual comparison of forecasts smpl @all graph grfore.line total_cases Gompertz_curve_f Gompertz_obs_f Gompertz_obs_intervention_f Gompertz_obs_dynamic_f grfore.draw(shade, bottom, @rgb(255,235,204)) {%last} {%end} 'Compare the forecasts of Dynamic Gompertz model before and after the beginning of normalization period series total_cases_before = @recode(@date<@dateval("6/1/2020"),total_cases,NA) tsepigrowth(ser="total_cases_before", growth="1", model="3", sample="3/12/2020 12/31/2020") ssmodel03.makesignals(t=smooth) Gompertz_obs_dynamic_g_before smpl "5/31/2020" %end series Gompertz_obs_dynamic_f_before = @recode(total_cases_before<>NA,total_cases_before(-1)*(1+@exp(Gompertz_obs_dynamic_g_before)),Gompertz_obs_dynamic_f_before(-1)*(1+@exp(Gompertz_obs_dynamic_g_before))) smpl @all graph grfore2.line total_cases total_cases_before Gompertz_obs_dynamic_f_before Gompertz_obs_dynamic_f grfore2.draw(shade, bottom, @rgb(255,235,204)) 6/1/2020 {%last}