/* This program generates some Monte Carlo data for various permanent intervention models. */ options pagesize=60 linesize=74 nodate; goptions device=win gsfname=plot rotate=landscape gsfmode=append target=hpljs3; *border; /* This segment of the program generates the deterministic part of an abrupt start, permanent change intervention model. */ data mc1; do t = 0 to 100; S = (t ge 25); y = 4*S; if t > 0 then output; end; keep t y S; proc gplot data=mc1; symbol v=dot c=black i=join h=.8; title1 'Monte Carlo Intervention Data with Abrupt Start, Permanent Change'; title2 '(With No Noise)'; axis1 order=(0 to 100 by 10) label=(f=duplex 'Obs'); axis2 order=(-10 to 10 by 2) label=(f=duplex 'Intervention Series'); plot y*t / haxis=axis1 vaxis=axis2; run; /* This segment of the program adds noise to the previous deterministic intervention. */ data mc2; iseed = 13333; do t = 0 to 100; S = (t ge 25); a = rannor(iseed); y = 4*S + a; if t > 0 then output; end; keep t y S; proc gplot data=mc2; symbol v=dot c=black i=join h=.8; title1 'Monte Carlo Intervention Data with Abrupt Start, Permanent Change'; title2 '(With Noise)'; axis1 order=(0 to 100 by 10) label=(f=duplex 'Obs'); axis2 order=(-10 to 10 by 2) label=(f=duplex 'Intervention Series'); plot y*t / haxis=axis1 vaxis=axis2; run; /* Here we estimate the above model. */ proc arima data = mc2; identify var = y crosscorr = (S) noprint; estimate input = (S) plot; estimate p = 1 input = (S); estimate q = 1 input = (S); run; /* This segment of the program generates the deterministic part of an slow start, permanent change intervention model. */ data mc3; y1 = 0; do t = 0 to 100; S = (t ge 25); y2 = 0.9*y1 + 1*S; if t > 0 then output; y1 = y2; end; keep t y2 S; proc gplot data=mc3; symbol v=dot c=black i=join h=.8; title1 'Monte Carlo Intervention Data with Slow Start, Permanent Change'; title2 '(With No Noise)'; axis1 order=(0 to 100 by 10) label=(f=duplex 'Obs'); axis2 order=(-6 to 16 by 2) label=(f=duplex 'Intervention Series'); plot y2*t / haxis=axis1 vaxis=axis2; run; /* This segment of the program adds noise to the previous deterministic intervention. */ data mc4; y1 = 0; iseed = 123; do t = 0 to 100; S = (t ge 25); a = rannor(iseed); y2 = 0.9*y1 + 1*S + 0.5*a; if t > 0 then output; y1 = y2; end; keep t y2 S; proc gplot data=mc4; symbol v=dot c=black i=join h=.8; title1 'Monte Carlo Intervention Data with Slow Start, Permanent Change'; title2 '(With Noise)'; axis1 order=(0 to 100 by 10) label=(f=duplex 'Obs'); axis2 order=(-6 to 16 by 2) label=(f=duplex 'Intervention Series'); plot y2*t / haxis=axis1 vaxis=axis2; run; /* Here we estimate the above model. */ proc arima data = mc4; identify var = y2 crosscorr = (S) noprint; estimate input = (/(1)S) plot; estimate p = 1 input = (/(1)S); estimate q = 1 input = (/(1)S); /* These last 4 models are overfits of the second model. */ estimate p = 2 input = (/(1)S); estimate p = 1 q = 1 input = (/(1)S); estimate p = 1 input = (/(1,2)S); estimate p = 1 input = ((1)/(1)S); /* The second model is the preferred one. */ run; /* This segment of the program generates the deterministic part of an slow start, permanent change intervention model with "wave approach". */ data mc5; y1 = 0; y2 = 0; do t = 0 to 100; S = (t ge 25); y3 = 1.6*y2 - 0.8*y1 + 1*S; if t > 0 then output; y1 = y2; y2 = y3; end; keep t y3 S; proc gplot data=mc5; symbol v=dot c=black i=join h=.8; title1 'Monte Carlo Intervention Data with Slow Start, Permanent Change'; title2 'Wave Approach (With No Noise)'; axis1 order=(0 to 100 by 10) label=(f=duplex 'Obs'); axis2 order=(-6 to 14 by 2) label=(f=duplex 'Intervention Series'); plot y3*t / haxis=axis1 vaxis=axis2; run; /* This segment of the program adds noise to the previous deterministic intervention. */ data mc6; y1 = 0; y2 = 0; iseed = 19611; do t = 0 to 100; S = (t ge 25); a = rannor(iseed); y3 = 1.6*y2 - 0.8*y1 + 1*S + 0.05*a; if t > 0 then output; y1 = y2; y2 = y3; end; keep t y3 S; proc gplot data=mc6; symbol v=dot c=black i=join h=.8; title1 'Monte Carlo Intervention Data with Slow Start, Permanent Change'; title2 'Wave Approach (With Noise)'; axis1 order=(0 to 100 by 10) label=(f=duplex 'Obs'); axis2 order=(-6 to 14 by 2) label=(f=duplex 'Intervention Series'); plot y3*t / haxis=axis1 vaxis=axis2; run; /* Here we estimate the above model. */ proc arima data = mc6; identify var = y3 crosscorr = (S) noprint; estimate input = (/(1,2)S) plot; estimate p = 1 input = (/(1,2)S); estimate q = 1 input = (/(1,2)S); estimate p = 2 input = (/(1,2)S); /* The next 4 models are overfits of the fourth model. */ estimate p = 3 input = (/(1,2)S); estimate p = 2 q = 1 input = (/(1,2)S); estimate p = 2 input = (/(1,2,3)S); estimate p = 2 input = ((1)/(1,2)S); /* The fourth model is the preferred one. */ run;