forked from sassoftware/covid-19-sas
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathsimseir.sas
121 lines (96 loc) · 2.27 KB
/
simseir.sas
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
ods graphics on;
%let N=50e6;
%let tau=5.1;
%let Rho0=2.5;
%let sigma=.9;
%let start='1mar2020'd;
%let end='1sep2020'd;
data timepts;
s = &N; e = 0; v = 1; r = 0;
do date = &start to &end;
output;
s = .; e = .; v = .; r = .;
end;
run;
data s;
_name_ = "cases";
cases = 0;
run;
data r0cov;
length _type_ $ 8 _name_ $ 8;
_type_ = "OLS";
_name_ = "";
R0 = &Rho0;
output;
_name_ = "R0";
R0 = 1;
output;
run;
data r0tcov;
length _type_ $ 8 _name_ $ 8;
_type_ = "OLS";
_name_ = "";
R0 = &Rho0;
tau = τ
output;
_name_ = "R0";
R0 = 0.5; tau = 0;
output;
_name_ = "tau";
R0 = 0; tau = 2;
output;
run;
title 'SEIR Monte Carlo simulation';
proc tmodel data=timepts;
endo s &N e 0 v 1 r 0;
/* parameters of interest */
parms R0 &Rho0 tau τ
/* fixed values */
control N &N
sigma σ
/* coefficient parameterizations */
gamma = 1/tau;
beta = R0*gamma/N;
/* Differential equations */
dert.s = - beta*s*v;
dert.e = beta*s*v - sigma*e;
dert.v = sigma*e - gamma*v;
dert.r = gamma*v;
cases = v + r;
/* monte carlo simulation; r0 */
solve cases / time=date outpredict out=mcsimr0(rename=v=i) seed=1
random=5 sdata=s estdata=r0cov;
outvars R0 tau;
/* monte carlo simulation: r0 and tau */
solve cases / time=date outpredict out=mcsimr0t(rename=v=i) seed=1
random=10 sdata=s estdata=r0tcov;
outvars R0 tau;
quit;
/* Hundred thousands Fomat */
proc format;
picture hundthou 0-high = '000.0' (mult=0.00001);
quit;
data mcsimr0;
set mcsimr0;
format date date.;
run;
proc print data=mcsimr0(where=(date=&start));
format s e i r hundthou.;
run;
proc sgplot data=mcsimr0 noautolegend;
series x =date y=i / group=_rep_ lineattrs=(pattern=solid thickness=2);
yaxis label='Number of Cases (x100,000)' valuesformat=hundthou.;
xaxis label='Date';
run;
data mcsimr0t;
set mcsimr0t;
format date date.;
run;
proc print data=mcsimr0t(where=(date=&start));
format s e i r hundthou.;
run;
proc sgplot data=mcsimr0t noautolegend;
series x =date y=i / group=_rep_ lineattrs=(pattern=solid thickness=2);
yaxis label='Number of Cases (x100,000)' valuesformat=hundthou.;
xaxis label='Date';
run;