library(data.table)
library(forecast) #autoplot

Ontario <- fread("Cases of opioid-related morbidity and mortality, Ontario, 2003 - 01 – 2018 - 09.csv")

Toronto_Public_Health <- fread("Cases of opioid-related morbidity and mortality, Toronto Public Health, 2003 - 01 – 2018 - 09.csv")

str(Toronto_Public_Health)

Classes ‘data.table’ and 'data.frame':  189 obs. of  8 variables:
 $ Year - Month             : chr  "2003 - 01" "2003 - 02" "2003 - 03" "2003 - 04" ...
 $ Population               : int  2600854 2600228 2599602 2598977 2598351 2597725 2597099 2596801 2596503 2596205 ...
 $ Rate of ED visits        : num  1 0.5 0.8 0.8 1.1 0.8 1.1 0.9 1.2 1 ...
 $ Count of ED visits       : int  25 13 22 22 28 20 28 24 30 27 ...
 $ Rate of Hospitalizations : num  0.6 0.3 0.5 0.6 0.6 0.5 0.5 0.7 0.4 0.6 ...
 $ Count of Hospitalizations: int  15 9 12 15 15 14 12 18 11 16 ...
 $ Rate of Deaths           : chr  "NA" "NA" "NA" "NA" ...
 $ Count of Deaths          : chr  "NA" "NA" "NA" "NA" ...
 - attr(*, ".internal.selfref")=<externalptr> 

Toronto_Central_LHIN <- fread("Cases of opioid-related morbidity and mortality, Toronto Central LHIN, 2003 - 01 – 2018 - 09.csv")

Ontario_Count_ED_visits.ts <- ts(Ontario$`Count of ED visits`, start = 2003, frequency = 12)

Toronto_Public_Health_Count_ED_visits.ts <- ts(Toronto_Public_Health$`Count of ED visits`, start = 2003, frequency = 12)

Toronto_Central_LHIN_Count_ED_visits.ts <- ts(Toronto_Central_LHIN$`Count of ED visits`, start = 2003, frequency = 12)

Ontario_Count_ED_visits.ts.components <- decompose(Ontario_Count_ED_visits.ts)

Toronto_Public_Health_Count_ED_visits.ts.components <- decompose(Toronto_Public_Health_Count_ED_visits.ts)

plot(Ontario_Count_ED_visits.ts.components)

plot(Toronto_Public_Health_Count_ED_visits.ts.components)

Toronto_Central_LHIN_Count_ED_visits.ts.components <- decompose(Toronto_Central_LHIN_Count_ED_visits.ts)

plot(Toronto_Central_LHIN_Count_ED_visits.ts.components)

Toronto_Public_Health_Count_ED_visits.ts.arima <- auto.arima(Toronto_Public_Health_Count_ED_visits.ts, D=1)

Toronto_Public_Health_Count_ED_visits.ts.arima

Series: Toronto_Public_Health_Count_ED_visits.ts 
ARIMA(3,1,1)(1,1,1)[12] 

Coefficients:
         ar1     ar2      ar3      ma1     sar1     sma1
      0.7587  0.0642  -0.2713  -0.7898  -0.2953  -0.6304
s.e.  0.1014  0.0939   0.0771   0.0834   0.1300   0.1191

sigma^2 estimated as 177.3:  log likelihood=-708.23
AIC=1430.46   AICc=1431.13   BIC=1452.65

Toronto_Public_Health_Count_ED_visits.ts.arima.forecast <- forecast(Toronto_Public_Health_Count_ED_visits.ts.arima, level = c(95), h = 24)

autoplot(Toronto_Public_Health_Count_ED_visits.ts.arima.forecast, main = "Monthly Emergency Department (ED) Visits")

Toronto_Public_Health_Rate_of_ED_visits.ts <- ts(Toronto_Public_Health$`Rate of ED visits`, start = 2003, frequency = 12)

Toronto_Public_Health_Rate_of_ED_visits.ts.components <- decompose(Toronto_Public_Health_Rate_of_ED_visits.ts)

plot(Toronto_Public_Health_Rate_of_ED_visits.ts.components)

Toronto_Public_Health_Rate_of_ED_visits.ts.arima <- auto.arima(Toronto_Public_Health_Rate_of_ED_visits.ts, D=1)

Toronto_Public_Health_Rate_of_ED_visits.ts.arima

Series: Toronto_Public_Health_Rate_of_ED_visits.ts 
ARIMA(0,1,0)(1,1,1)[12] 

Coefficients:
         sar1     sma1
      -0.2606  -0.6268
s.e.   0.1348   0.1218

sigma^2 estimated as 0.2475:  log likelihood=-131.09
AIC=268.18   AICc=268.32   BIC=277.7

Toronto_Public_Health_Rate_of_ED_visits.ts.arima.forecast <- forecast(Toronto_Public_Health_Rate_of_ED_visits.ts.arima, level = c(95), h = 24)

autoplot(Toronto_Public_Health_Rate_of_ED_visits.ts.arima.forecast, main = "Monthly Emergency Department (ED) Visits")

Ontario_Count_ED_visits.ts.arima <- auto.arima(Ontario_Count_ED_visits.ts, D=1)

Ontario_Count_ED_visits.ts.arima

Series: Ontario_Count_ED_visits.ts 
ARIMA(0,1,0)(0,1,2)[12] 

Coefficients:
         sma1    sma2
      -0.8092  0.1603
s.e.   0.0885  0.0816

sigma^2 estimated as 1473:  log likelihood=-894.99
AIC=1795.98   AICc=1796.12   BIC=1805.49

Ontario_Count_ED_visits.ts.arima.forecast <- forecast(Ontario_Count_ED_visits.ts.arima, level = c(95), h = 24)

autoplot(Ontario_Count_ED_visits.ts.arima.forecast, main = "Monthly Emergency Department (ED) Visits")

LS0tDQp0aXRsZTogIlRvcm9udG8gYW5kIE9udGFyaW8gT3Bpb2lkIE1vcmJpZGl0eSBUaW1lc2VyaWVzIg0Kb3V0cHV0OiBodG1sX25vdGVib29rDQotLS0NCg0KIyBBdXRob3INCg0KKiBKb3JkYW4gQmVsbA0KKiBKdWx5IDExLCAyMDE5DQoqIDxodHRwczovL2pvcmRhbmJlbGwyMzU3LmdpdGh1Yi5pby9PcGlvaWRfVGltZXNlcmllcy5uYi5odG1sPg0KDQojIFB1YmxpYyBIZWFsdGggT250YXJpbw0KDQpQdWJsaWMgSGVhbHRoIE9udGFyaW8sIFtJbnRlcmFjdGl2ZSBPcGlvaWQgVG9vbF0oaHR0cHM6Ly93d3cucHVibGljaGVhbHRob250YXJpby5jYS9lbi9kYXRhLWFuZC1hbmFseXNpcy9zdWJzdGFuY2UtdXNlL2ludGVyYWN0aXZlLW9waW9pZC10b29sIy9hZ2VTZXgpDQoNCk9waW9pZC1yZWxhdGVkIG1vcmJpZGl0eSBhbmQgbW9ydGFsaXR5IGluIE9udGFyaW8NCg0KYGBge3J9DQpsaWJyYXJ5KGRhdGEudGFibGUpDQpsaWJyYXJ5KGZvcmVjYXN0KSAjYXV0b3Bsb3QNCmBgYA0KDQpgYGB7cn0NCk9udGFyaW8gPC0gZnJlYWQoIkNhc2VzIG9mIG9waW9pZC1yZWxhdGVkIG1vcmJpZGl0eSBhbmQgbW9ydGFsaXR5LCBPbnRhcmlvLCAyMDAzIC0gMDEg4oCTIDIwMTggLSAwOS5jc3YiKQ0KYGBgDQoNCmBgYHtyfQ0KVG9yb250b19QdWJsaWNfSGVhbHRoIDwtIGZyZWFkKCJDYXNlcyBvZiBvcGlvaWQtcmVsYXRlZCBtb3JiaWRpdHkgYW5kIG1vcnRhbGl0eSwgVG9yb250byBQdWJsaWMgSGVhbHRoLCAyMDAzIC0gMDEg4oCTIDIwMTggLSAwOS5jc3YiKQ0KYGBgDQoNCmBgYHtyfQ0Kc3RyKFRvcm9udG9fUHVibGljX0hlYWx0aCkNCmBgYA0KDQpgYGB7cn0NClRvcm9udG9fQ2VudHJhbF9MSElOIDwtIGZyZWFkKCJDYXNlcyBvZiBvcGlvaWQtcmVsYXRlZCBtb3JiaWRpdHkgYW5kIG1vcnRhbGl0eSwgVG9yb250byBDZW50cmFsIExISU4sIDIwMDMgLSAwMSDigJMgMjAxOCAtIDA5LmNzdiIpDQpgYGANCg0KIyBUcmVuZCBhbmQgU2Vhc29uYWxpdHkgRGVjb21wb3NpdGlvbg0KDQpgYGB7cn0NCk9udGFyaW9fQ291bnRfRURfdmlzaXRzLnRzIDwtIHRzKE9udGFyaW8kYENvdW50IG9mIEVEIHZpc2l0c2AsIHN0YXJ0ID0gMjAwMywgZnJlcXVlbmN5ID0gMTIpDQpgYGANCg0KYGBge3J9DQpUb3JvbnRvX1B1YmxpY19IZWFsdGhfQ291bnRfRURfdmlzaXRzLnRzIDwtIHRzKFRvcm9udG9fUHVibGljX0hlYWx0aCRgQ291bnQgb2YgRUQgdmlzaXRzYCwgc3RhcnQgPSAyMDAzLCBmcmVxdWVuY3kgPSAxMikNCmBgYA0KDQpgYGB7cn0NClRvcm9udG9fQ2VudHJhbF9MSElOX0NvdW50X0VEX3Zpc2l0cy50cyA8LSB0cyhUb3JvbnRvX0NlbnRyYWxfTEhJTiRgQ291bnQgb2YgRUQgdmlzaXRzYCwgc3RhcnQgPSAyMDAzLCBmcmVxdWVuY3kgPSAxMikNCmBgYA0KDQpgYGB7cn0NCk9udGFyaW9fQ291bnRfRURfdmlzaXRzLnRzLmNvbXBvbmVudHMgPC0gZGVjb21wb3NlKE9udGFyaW9fQ291bnRfRURfdmlzaXRzLnRzKQ0KYGBgDQoNCmBgYHtyfQ0KVG9yb250b19QdWJsaWNfSGVhbHRoX0NvdW50X0VEX3Zpc2l0cy50cy5jb21wb25lbnRzIDwtIGRlY29tcG9zZShUb3JvbnRvX1B1YmxpY19IZWFsdGhfQ291bnRfRURfdmlzaXRzLnRzKQ0KYGBgDQoNCmBgYHtyfQ0KcGxvdChPbnRhcmlvX0NvdW50X0VEX3Zpc2l0cy50cy5jb21wb25lbnRzKQ0KYGBgDQoNCg0KYGBge3J9DQpwbG90KFRvcm9udG9fUHVibGljX0hlYWx0aF9Db3VudF9FRF92aXNpdHMudHMuY29tcG9uZW50cykNCmBgYA0KDQpgYGB7cn0NClRvcm9udG9fQ2VudHJhbF9MSElOX0NvdW50X0VEX3Zpc2l0cy50cy5jb21wb25lbnRzIDwtIGRlY29tcG9zZShUb3JvbnRvX0NlbnRyYWxfTEhJTl9Db3VudF9FRF92aXNpdHMudHMpDQpgYGANCg0KYGBge3J9DQpwbG90KFRvcm9udG9fQ2VudHJhbF9MSElOX0NvdW50X0VEX3Zpc2l0cy50cy5jb21wb25lbnRzKQ0KYGBgDQoNCiMgU0FSSU1BIG1vZGVscyBhbmQgZm9yZWNhc3RzDQoNCkQ9MSBmb3JjZXMgc2Vhc29uYWxpdHkNCg0KYGBge3J9DQpUb3JvbnRvX1B1YmxpY19IZWFsdGhfQ291bnRfRURfdmlzaXRzLnRzLmFyaW1hIDwtIGF1dG8uYXJpbWEoVG9yb250b19QdWJsaWNfSGVhbHRoX0NvdW50X0VEX3Zpc2l0cy50cywgRD0xKQ0KDQpUb3JvbnRvX1B1YmxpY19IZWFsdGhfQ291bnRfRURfdmlzaXRzLnRzLmFyaW1hDQpgYGANCg0KYGBge3J9DQpUb3JvbnRvX1B1YmxpY19IZWFsdGhfQ291bnRfRURfdmlzaXRzLnRzLmFyaW1hLmZvcmVjYXN0IDwtIGZvcmVjYXN0KFRvcm9udG9fUHVibGljX0hlYWx0aF9Db3VudF9FRF92aXNpdHMudHMuYXJpbWEsIGxldmVsID0gYyg5NSksIGggPSAyNCkNCg0KYXV0b3Bsb3QoVG9yb250b19QdWJsaWNfSGVhbHRoX0NvdW50X0VEX3Zpc2l0cy50cy5hcmltYS5mb3JlY2FzdCwgbWFpbiA9ICJNb250aGx5IEVtZXJnZW5jeSBEZXBhcnRtZW50IChFRCkgVmlzaXRzIikNCmBgYA0KDQpgYGB7cn0NClRvcm9udG9fUHVibGljX0hlYWx0aF9SYXRlX29mX0VEX3Zpc2l0cy50cyA8LSB0cyhUb3JvbnRvX1B1YmxpY19IZWFsdGgkYFJhdGUgb2YgRUQgdmlzaXRzYCwgc3RhcnQgPSAyMDAzLCBmcmVxdWVuY3kgPSAxMikNCmBgYA0KDQpgYGB7cn0NClRvcm9udG9fUHVibGljX0hlYWx0aF9SYXRlX29mX0VEX3Zpc2l0cy50cy5jb21wb25lbnRzIDwtIGRlY29tcG9zZShUb3JvbnRvX1B1YmxpY19IZWFsdGhfUmF0ZV9vZl9FRF92aXNpdHMudHMpDQoNCnBsb3QoVG9yb250b19QdWJsaWNfSGVhbHRoX1JhdGVfb2ZfRURfdmlzaXRzLnRzLmNvbXBvbmVudHMpDQpgYGANCg0KYGBge3J9DQpUb3JvbnRvX1B1YmxpY19IZWFsdGhfUmF0ZV9vZl9FRF92aXNpdHMudHMuYXJpbWEgPC0gYXV0by5hcmltYShUb3JvbnRvX1B1YmxpY19IZWFsdGhfUmF0ZV9vZl9FRF92aXNpdHMudHMsIEQ9MSkNCg0KVG9yb250b19QdWJsaWNfSGVhbHRoX1JhdGVfb2ZfRURfdmlzaXRzLnRzLmFyaW1hDQpgYGANCg0KYGBge3J9DQpUb3JvbnRvX1B1YmxpY19IZWFsdGhfUmF0ZV9vZl9FRF92aXNpdHMudHMuYXJpbWEuZm9yZWNhc3QgPC0gZm9yZWNhc3QoVG9yb250b19QdWJsaWNfSGVhbHRoX1JhdGVfb2ZfRURfdmlzaXRzLnRzLmFyaW1hLCBsZXZlbCA9IGMoOTUpLCBoID0gMjQpDQoNCmF1dG9wbG90KFRvcm9udG9fUHVibGljX0hlYWx0aF9SYXRlX29mX0VEX3Zpc2l0cy50cy5hcmltYS5mb3JlY2FzdCwgbWFpbiA9ICJNb250aGx5IEVtZXJnZW5jeSBEZXBhcnRtZW50IChFRCkgVmlzaXRzIikNCmBgYA0KDQpgYGB7cn0NCk9udGFyaW9fQ291bnRfRURfdmlzaXRzLnRzLmFyaW1hIDwtIGF1dG8uYXJpbWEoT250YXJpb19Db3VudF9FRF92aXNpdHMudHMsIEQ9MSkNCg0KT250YXJpb19Db3VudF9FRF92aXNpdHMudHMuYXJpbWENCmBgYA0KDQpgYGB7cn0NCk9udGFyaW9fQ291bnRfRURfdmlzaXRzLnRzLmFyaW1hLmZvcmVjYXN0IDwtIGZvcmVjYXN0KE9udGFyaW9fQ291bnRfRURfdmlzaXRzLnRzLmFyaW1hLCBsZXZlbCA9IGMoOTUpLCBoID0gMjQpDQoNCmF1dG9wbG90KE9udGFyaW9fQ291bnRfRURfdmlzaXRzLnRzLmFyaW1hLmZvcmVjYXN0LCBtYWluID0gIk1vbnRobHkgRW1lcmdlbmN5IERlcGFydG1lbnQgKEVEKSBWaXNpdHMiKQ0KYGBgDQoNCg==