Science, Irving K. Barber Faculty of (Okanagan)
Computer Science, Mathematics, Physics and Statistics, Department of (Okanagan)
DSpace
UBCO
Snyman, Simon
Master of Science - MSc
University of British Columbia
Exponential smoothing methods offer several tools for forecasting and simulating varying time series patterns, including historical wildfire counts. The Box-Cox transform, ARMA errors, trend and seasonal components (BATS), and the trigonometric BATS model (TBATS) are exponential smoothing techniques capable of modeling complex seasonality, non-integer frequencies, and more. The BATS and TBATS frameworks support short- and intermediate-term forecasting and simulation, while also providing a foundation for model-based time series clustering of wildfire counts. The application in this thesis of the BATS and TBATS models aims to support fire agencies in forecasting and simulating wildfire occurrences across Canada on an increased scale that provides alternative methods for future use. Variable length Markov chains (VLMCs) are an extension of traditional fixed-order Markov chains that can reduce model complexity and improve computational efficiency. To assess the practicality and effectiveness of VLMCs in modeling wildfire causes in Canada, they are evaluated in terms of forecasting on both balanced and unbalanced data. Furthermore, VLMCs are used to model the cluster assignment results on the weekly wildfire count time series data over a range of years. We provide preliminary information on clustering weekly wildfire counts using a model-based approach built from TBATS models that can help identify broad characteristics for an upcoming fire season through extensions of Markov chains.
Clustering Wildfire Occurrence TimeSeriesbySimon Marcus SnymanB.Sc. Mathematics, The University of British Columbia, 2023A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFMASTER OF SCIENCEinTHE COLLEGE OF GRADUATE STUDIES(Mathematics)THE UNIVERSITY OF BRITISH COLUMBIA(Okanagan)August 2025© Simon Marcus Snyman, 2025The following individuals certify that they have read, and recommendto the College of Graduate Studies for acceptance, a thesis/dissertation en-titled:Clustering Wildfire Occurrence Time Seriessubmitted by Simon Marcus Snyman in partial fulfilment of the require-ments of the degree of Master of ScienceDr. W. John Braun, Irving K. Barber Faculty of ScienceCo-supervisorDr. Lengyi Han, Irving K. Barber Faculty of ScienceCo-supervisorDr. Irene Vrbik, Irving K. Barber Faculty of ScienceSupervisory Committee MemberDr. Michael Noonan, Irving K. Barber Faculty of ScienceSupervisory Committee MemberDr. Frederic Menard, Irving K. Barber Faculty of ScienceUniversity ExamineriiAbstractExponential smoothing methods offer several tools for forecasting andsimulating varying time series patterns, including historical wildfire counts.The Box-Cox transform, ARMA errors, trend and seasonal components(BATS), and the trigonometric BATS model (TBATS) are exponential smooth-ing techniques capable of modeling complex seasonality, non-integer frequen-cies, and more. The BATS and TBATS frameworks support short- andintermediate-term forecasting and simulation, while also providing a foun-dation for model-based time series clustering of wildfire counts.The application in this thesis of the BATS and TBATS models aimsto support fire agencies in forecasting and simulating wildfire occurrencesacross Canada on an increased scale that provides alternative methods forfuture use.Variable length Markov chains (VLMCs) are an extension of traditionalfixed-order Markov chains that can reduce model complexity and improvecomputational efficiency. To assess the practicality and effectiveness ofVLMCs in modeling wildfire causes in Canada, they are evaluated in termsof forecasting on both balanced and unbalanced data. Furthermore, VLMCsare used to model the cluster assignment results on the weekly wildfire counttime series data over a range of years. We provide preliminary informationon clustering weekly wildfire counts using a model-based approach built fromTBATS models that can help identify broad characteristics for an upcomingfire season through extensions of Markov chains.iiiLay SummaryTime series data are ordered observations collected over time and arecommonly used in areas such as finance, economics, and environmental sci-ence. One powerful use of time series data is to forecast and simulate futureobservations, which can provide valuable insight for real-world problems.Clustering is a method used to group data based on similarities. In thecontext of time series, clustering helps classify similar patterns based oncharacteristics such as shape or statistical features. This allows researchersto discover relationships and trends that might not be immediately obvious.This thesis focuses on forecasting, simulating, and clustering wildfirecounts and categorical wildfire causes throughout Canada. The number ofwildfires over time can be represented as a time series, while the causes ofwildfires can be treated as a time series of categorical outcomes. This frame-work offers new opportunities to gain insight into the behavior of wildfiresthroughout Canada.To apply time series forecasting, simulation, and clustering to wild-fire count data, exponential smoothing methods are explored and appliedacross several Canadian regions. The results demonstrate that these meth-ods can handle complex seasonal patterns and support accurate short- tointermediate-term forecasting. They also provide a foundation for a model-based clustering approach for time series of wildfires. In addition, wildfirecause data are analyzed using an extension of Markov chains known as vari-able length Markov chains (VLMCs). These analyses of wildfire causes offera more efficient approach for forecasting categorical wildfire outcomes thantraditional Markov models. The resulting model can be used to identifycharacteristics of the most likely upcoming fire season, given historical data.ivPrefaceThis thesis was written by Simon Snyman under the co-supervision ofDr. W. John Braun and Dr. Lengyi Han. The research problems addressedemerged through various conversations with Dr. Braun and Dr. Han, bothof whom provided valuable guidance and insightful feedback throughout theprocess.Simon Snyman was responsible for initiating, developing, and conduct-ing the data analysis, forecasting, simulations, and clustering procedurespresented in this work. Generative AI was not used for this work. Themanuscript was constructed by Simon Snyman through a combination oforiginal work for this thesis and research conferences.The wildfire count simulation, forecasting and clustering research wascarried out over a period of 12 months, in collaboration with Dr. Braun andDr. Han. Preliminary results were presented by Simon Snyman at the 2024meeting of the International Society of Environmetrics in Adelaide, Australia(9 December 2024) and at the Statistical Society of Canada meeting inSaskatoon, Saskatchewan (27 May 2025).The wildfire causes analysis was performed independently by Simon Sny-man, with feedback provided by Dr. Braun and Dr. Eric Foxall. This compo-nent of the research originated as a graduate course project at the Universityof British Columbia Okanagan.vTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiLay Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . viList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ixList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . xiiiDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xivChapter 1: Introduction . . . . . . . . . . . . . . . . . . . . . . . 11.1 Time Series . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.2 Forecasting and Clustering Wildfire Time Series . . . . . . . . 2Chapter 2: Background . . . . . . . . . . . . . . . . . . . . . . . 42.1 Time Series Foundations . . . . . . . . . . . . . . . . . . . . . 42.1.1 Time Series Properties . . . . . . . . . . . . . . . . . . 42.1.2 Autoregressive Models . . . . . . . . . . . . . . . . . . 52.1.3 Simple Exponential Smoothing Methods . . . . . . . . 62.1.4 Modified Models: TBATS and BATS Models . . . . . 92.2 State Space Models . . . . . . . . . . . . . . . . . . . . . . . . 102.2.1 Linear State Space Models . . . . . . . . . . . . . . . 112.2.2 Linear Gaussian State Space Models . . . . . . . . . . 122.2.3 State Estimation Techniques . . . . . . . . . . . . . . 132.2.4 State Space Innovations: BATS and TBATS . . . . . 14viTABLE OF CONTENTS2.3 Clustering Time Series Models . . . . . . . . . . . . . . . . . 162.3.1 Mixtures of Linear Gaussian State Space Models . . . 162.4 Markov Chains . . . . . . . . . . . . . . . . . . . . . . . . . . 182.4.1 Variable Length Markov Chains . . . . . . . . . . . . . 192.4.2 Variable Length Memory, Context Functions, and Con-text Trees . . . . . . . . . . . . . . . . . . . . . . . . . 202.5 Wildfire Occurrences in Canada . . . . . . . . . . . . . . . . . 222.5.1 Wildfire Occurrences Dataset . . . . . . . . . . . . . . 23Chapter 3: Methodology . . . . . . . . . . . . . . . . . . . . . . 273.1 Wildfire Occurrence Modeling . . . . . . . . . . . . . . . . . . 273.1.1 Dataset Description and Challenges . . . . . . . . . . 273.1.2 Model Selection and Fitting . . . . . . . . . . . . . . . 283.1.3 Residual Diagnostics: TBATS . . . . . . . . . . . . . . 293.1.4 Forecasting and Simulating Wildfire Time Series . . . 303.2 Clustering Wildfire Time Series . . . . . . . . . . . . . . . . . 313.2.1 Clustering Method and Feature Extraction . . . . . . 313.2.2 Mclust and the EM Algorithm . . . . . . . . . . . . . 333.2.3 Fitting VLMCs . . . . . . . . . . . . . . . . . . . . . . 35Chapter 4: Applications . . . . . . . . . . . . . . . . . . . . . . . 374.1 Simulation Studies . . . . . . . . . . . . . . . . . . . . . . . . 374.1.1 TBATS Simulations for Wildfire Data . . . . . . . . . 374.1.2 VLMC Wildfire Occurrence Data . . . . . . . . . . . . 414.2 Real-World Applications . . . . . . . . . . . . . . . . . . . . . 444.2.1 Predicting Wildfire Causes using VLMCs . . . . . . . 444.2.2 VLMC versus First Order Markov Chain: WildfireCause Class Imbalance . . . . . . . . . . . . . . . . . . 464.2.3 Intermediate-Term Forecasting of Wildfires . . . . . . 484.2.4 Clustering Wildfire Time Series . . . . . . . . . . . . . 50Chapter 5: Conclusion . . . . . . . . . . . . . . . . . . . . . . . . 565.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 565.2 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . 57Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58Appendices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62Appendix A: Tables . . . . . . . . . . . . . . . . . . . . . . . . . . 62Appendix B: Figures . . . . . . . . . . . . . . . . . . . . . . . . . . 69viiTABLE OF CONTENTSAppendix C: AR(1) Model Derivations . . . . . . . . . . . . . . . . 77Appendix D: BATS and TBATS Model Derivations . . . . . . . . 79viiiList of TablesTable 4.1 VLMC model summary of 500 simulated cluster as-signments. . . . . . . . . . . . . . . . . . . . . . . . . . 42Table 4.2 Comparison of VLMC first order Markov chain to anindependent sampling baseline by region for wildfirecause in Alberta. Metrics include out-of-sample accu-racy and macro-averaged F1 score. . . . . . . . . . . . 48Table 4.3 Cluster assignments for BC wildfire occurrence timeseries from 1950-2023. . . . . . . . . . . . . . . . . . . 54Table 5.1 Notable variables and descriptions from the CanadianNational Fire Database (CNFDB) dataset. . . . . . . 62Table 5.2 Components of BATS and TBATS model fits. TBATSincludes additional components for the seasonal com-ponents. . . . . . . . . . . . . . . . . . . . . . . . . . . 63Table 5.3 Model summary of fitted BATS and TBATS modelsfor weekly wildfire counts in British Columbia (2018-2023). . . . . . . . . . . . . . . . . . . . . . . . . . . . 64Table 5.4 Comparison of BATS and TBATS out-of-sample modelperformance for weekly wildfire counts in British Columbia(2018-2023). Forecast accuracy is measured using MSEand MAE for each year in the test set. . . . . . . . . . 64Table 5.5 Comparison of BATS and TBATS out-of-sample modelperformance for weekly wildfire counts in Alberta (2018-2023). Forecast accuracy is measured using MSE andMAE for each year in the test set. . . . . . . . . . . . 65Table 5.6 Comparison of BATS and TBATS out-of-sample modelperformance for weekly wildfire counts in Saskatchewan(2018-2023). Forecast accuracy is measured using MSEand MAE for each year in the test set. . . . . . . . . . 65ixLIST OF TABLESTable 5.7 Out-of-sample VLMC model metrics for southern Al-berta by cause. Wildfire cause distribution for theSouth region: H: 15,957; N: 6,162. . . . . . . . . . . . 66Table 5.8 Out-of-sample VLMC model metrics for northwest Al-berta by cause. Wildfire cause distribution for theNorthwest region: H: 10,562; N: 10,298. . . . . . . . . 67Table 5.9 Out-of-sample VLMC model metrics for northeast Al-berta by cause. Wildfire cause distribution for theNortheast region: H: 4,031; N: 7,585. . . . . . . . . . . 68xList of FiguresFigure 1.1 Model-based clustering with varying time series pat-terns modeled by exponential smoothing methods. . . 3Figure 2.1 Context tree τc for X = {H,N,U}. . . . . . . . . . . 22Figure 2.2 Weekly wildfire occurrences in British Columbia in2003. . . . . . . . . . . . . . . . . . . . . . . . . . . . 25Figure 2.3 Weekly wildfire occurrences in British Columbia from1950-2023. . . . . . . . . . . . . . . . . . . . . . . . . 26Figure 3.1 Model-based clustering with varying time series pat-terns modeled by TBATS models. . . . . . . . . . . . 33Figure 4.1 TBATS fit and simulation of weekly wildfires in BritishColumbia, 2003. . . . . . . . . . . . . . . . . . . . . . 38Figure 4.2 TBATS simulation of weekly wildfires in British Columbia,2018–2023. . . . . . . . . . . . . . . . . . . . . . . . . 39Figure 4.3 TBATS simulation using 5 iterations with 95% CI forwildfires in British Columbia, 2018–2023. . . . . . . . 40Figure 4.4 TBATS simulation using 50 iterations with 95% CIfor wildfires in British Columbia, 2018–2023. . . . . . 40Figure 4.5 TBATS simulation using 500 iterations with 95% CIfor wildfires in British Columbia, 2018–2023. . . . . . 41Figure 4.6 Comparison of yearly BATS and TBATS forecast ac-curacy in BC from 2018-2023 using MSE and MAE. . 49Figure 4.7 Comparison of yearly BATS and TBATS forecast ac-curacy in Alberta from 2018-2023 using MSE and MAE. 50Figure 4.8 Comparison of yearly BATS and TBATS forecast ac-curacy in Saskatchewan from 2018-2023 using MSEand MAE.. . . . . . . . . . . . . . . . . . . . . . . . . 50Figure 4.9 TBATS simulated wildfire counts for Cluster 1 inBritish Columbia. . . . . . . . . . . . . . . . . . . . . 52xiLIST OF FIGURESFigure 4.10 TBATS simulated wildfire counts for Cluster 2 inBritish Columbia. . . . . . . . . . . . . . . . . . . . . 52Figure 4.11 TBATS simulated wildfire counts for Cluster 3 inBritish Columbia. . . . . . . . . . . . . . . . . . . . . 53Figure 5.1 TBATS simulation with 5 iterations: British Columbia,2018. . . . . . . . . . . . . . . . . . . . . . . . . . . . 69Figure 5.2 TBATS simulation with 5 iterations: British Columbia,2019. . . . . . . . . . . . . . . . . . . . . . . . . . . . 70Figure 5.3 TBATS simulation with 5 iterations: British Columbia,2020. . . . . . . . . . . . . . . . . . . . . . . . . . . . 70Figure 5.4 TBATS simulation with 5 iterations: British Columbia,2021. . . . . . . . . . . . . . . . . . . . . . . . . . . . 71Figure 5.5 TBATS simulation with 5 iterations: British Columbia,2022. . . . . . . . . . . . . . . . . . . . . . . . . . . . 71Figure 5.6 TBATS simulation with 5 iterations: British Columbia,2023. . . . . . . . . . . . . . . . . . . . . . . . . . . . 72Figure 5.7 TBATS simulation with 50 iterations: British Columbia,2018. . . . . . . . . . . . . . . . . . . . . . . . . . . . 72Figure 5.8 TBATS simulation with 50 iterations: British Columbia,2019. . . . . . . . . . . . . . . . . . . . . . . . . . . . 73Figure 5.9 TBATS simulation with 50 iterations: British Columbia,2020. . . . . . . . . . . . . . . . . . . . . . . . . . . . 73Figure 5.10 TBATS simulation with 50 iterations: British Columbia,2021. . . . . . . . . . . . . . . . . . . . . . . . . . . . 74Figure 5.11 TBATS simulation with 50 iterations: British Columbia,2022. . . . . . . . . . . . . . . . . . . . . . . . . . . . 74Figure 5.12 TBATS simulation with 50 iterations: British Columbia,2023. . . . . . . . . . . . . . . . . . . . . . . . . . . . 75Figure 5.13 Weekly wildfire occurrences in British Columbia from1998-2024. . . . . . . . . . . . . . . . . . . . . . . . . 75Figure 5.14 Wildfire Causes in Alberta (1959–2023). . . . . . . . . 76xiiAcknowledgementsI would like to thank my Mom and Dad for supporting me and believingin me over the years. All the conversations, family gatherings you organized,and golf kept me in good health during my studies. You pushed me to alwaysput forth my best effort and step out of my comfort zone, none of this wouldbe possible without you.Many thanks to my remaining family for their unconditional supportand encouragement throughout my academic journey so far. From Oumaand Oupa to my siblings, I thank you all for your unconditional supportover the years.I also would like to thank all my friends who have been there for mefrom the beginning of my undergraduate studies to this point. To my bestfriends Josh, Quinn, Kade, and Connor, I would like to thank you all forbeing great friends, supporting me, and providing me with friendships thatwill last a lifetime.To my girlfriend Ashley, I want to thank you for your unconditionalsupport and believing in me even when I did not. Through all of this, youmade me a better person both academically and professionally and I couldnot have finished any of this work without you.This research was funded in part by Wired NSF-NSERC Global Centreand a MITACS internship with FortisBC. I am grateful for the funding tosupport this research during my graduate studies.Finally, I would like to thank both my awesome co-supervisors Dr. W.John Braun and Dr. Lengyi Han for their guidance, support, and constantbelief. Throughout my master’s degree, Dr. Braun encouraged me to trynew ideas, topics, and instilled curiosity and interest in academics in methat I will carry forever. I also would like to thank Dr. Han for encouragingme during my undergraduate studies to pursue this work and believing inme. They both had an immense impact on my life and I would not be herewithout them.xiiiDedicationTo Ouma DalenexivChapter 1Introduction1.1 Time SeriesData points recorded over a time interval that are ordered are calledtime series [34]. The importance of sequential ordering of time series anal-ysis differentiates it from other known statistical methods [2]. Specifically,time series observations may depend on previous observations and otherspecific indexes in the sequence of observations, in contrast to independentobservations [2].There exist numerous possible patterns in time series data, such as whitenoise, cyclical trends, and seasonal patterns. Analysis of such patterns,enable effective modeling in fields such as economics, finance, engineering,and more. Seasonality in time series occurs when regular and predictablechanges are influenced by seasonal factors, such as the day of the weekor the month of the year [16]. Some time series with complex seasonalpatterns, including multiple seasonal periods, non-integer seasonality, andhigh-frequency seasonality, are common in real-world settings. The periodof a seasonal time series indicates when seasonal behavior tends to repeat,measured in units such as seconds, minutes, or hours.Topics to consider for different time series patterns include forecasting,simulating, clustering, and more. Specifically, various known exponentialsmoothing methods discussed in [8] will be applied in this thesis for differenttime series patterns and applications. Additionally, time series data can beexpressed as an ordered sequence of categorical observations over differenttime periods and regions. This may lead to forecasts or simulations that canbe useful in inferring and predicting future outcomes for various applications.The categorical time series methods applied in this thesis are inspired by[24].11.2. FORECASTING AND CLUSTERING WILDFIRE TIME SERIES1.2 Forecasting and Clustering Wildfire TimeSeriesThe focus of time series applications in this thesis will be on wildfiredata collected in Canada from the Fire Point and Polygon dataset from theCanadian National Fire Database (CNFDB) [6]. Specifically, the provincesof Alberta, British Columbia, and Saskatchewan will be evaluated in thiswork through modeling and clustering procedures. The main application ofexponential smoothing methods for time series data will be on the occurrenceof wildfires in western Canada. This thesis considers two complex exponen-tial smoothing methods from [8], where the effectiveness of modeling wildfireoccurrences in Canada is of interest to ultimately extend to a model-basedclustering approach. This is important as the number of wildfires in Canadais expected to increase in coming years [41]; Having a method that can helpidentify an upcoming class of wildfire occurrences over a calendar year canhelp with issues such as resource allocation, damage mitigation, and so on.Along with the model-based clustering approach for weekly wildfire occur-rences, we can identify the most likely class for an upcoming year. Thiscan provide information for fire agencies for upcoming fire seasons. Thesepreliminary model-based clustering results of the occurrence of wildfires inwestern Canada are presented with the aim of identifying which class afuture fire season could belong to. We employ exponential smoothing meth-ods from [8] and extensions of Markov chains in [24]. Figure 1.1 provides aschematic illustration of the model-based clustering approach in this thesiswith varying time series patterns.Model-based clustering results provide cluster assignments for weeklywildfire occurrences over a calendar year, which can be modeled with Markovchain methods. Our approach will be via the variable length Markov chainmethodology developed in [24]. Using the fitted Markov chain, informationfrom previous years can be used to match the class for the upcoming wild-fire season and provide inference on the most likely class for the next year.Given the class to which a future could belong, the model-based cluster-ing approach allows simulations and forecasts that can provide insight intofuture trends and patterns.21.2. FORECASTING AND CLUSTERING WILDFIRE TIME SERIESFigure 1.1: Model-based clustering illustration with varying time seriespatterns modelled by exponential smoothing methods on each cluster rep-resentation.3Chapter 2Background2.1 Time Series FoundationsTime series data are an ordered sequence of random variables that aresampled at equally spaced discrete points in a certain time interval. The in-spiration for these models came from applications in fields such as economics,business, and more [32]. Specifically, purposes such as the investigation ofthe dynamic relationship between variables, producing point or interval fore-casts, and understanding the dynamic structure of a process are all commonreal-world goals for time series analysis [32].The importance of time series applications is not limited to economics;specifically, social scientists, epidemiologists, and engineers should be inter-ested in time series models for investigating a range of problems [29]. Forexample, an epidemiologist could be concerned with the number of docu-mented COVID-19 cases in Canada during a specific time interval.The variety of applications and literature available for time series analysiscan be traced back to the early 1900s [32]. It provides a framework thathas motivated many fields to implement time series models to investigateproblems in both historical and modern contexts. This work will explorea wildfire application with different time series models and formulations.Specifically, the use of modelling inference and clustering will be the mainfocus of these applications in wildfire science.2.1.1 Time Series PropertiesIn this work, a time series will be denoted by {y1, y2, . . . , yn}, where tindicates the index. A typical way to express a time series isyt = f(Ψt−1) + et, (2.1)where Ψt−1 is in the σ−field that is generated by all previous values yn, and{et} is a sequence of independent identically distributed random variableswith zero mean and constant variance [32]. Historically, time series analysisfocuses on the function f(Ψt−1) and different errors {et} [32]. Moreover,42.1. TIME SERIES FOUNDATIONSdifferent forms of f and errors {et} can lead to different time series modelsto handle different types of data.A time series yt, as defined in Equation (2.1), possesses several key prop-erties, including mean function, variance function, autocovariance, and au-tocorrelation [7]. These properties at time t are defined as follows:µt = E[yt] (Mean function) (2.2)σ2t = Var(yt) (Variance function) (2.3)γk = Cov(yt, yt−k) (Autocovariance at lag k) (2.4)ρk =Cov(yt, yt−k)√Var(yt)Var(yt−k)(Autocorrelation at lag k). (2.5)The mean and variance of a time series yt is important to define stationarity;A time series is said to be (weakly) stationary if its mean and variance areconstant over time and the autocovariance γk depends only on k, not ontime t [32].Autocovariance measures the strength of the linear relationship betweenobservations at time t and t−k, while autocorrelation indicates the strengthof the correlation between these lags. Moreover, autocovariance and auto-correlation at lag k provide information on the linear dependence structureof a time series. Specifically, when the sequence of observations is indepen-dent, the autocorrelations at all lags are zero [13].2.1.2 Autoregressive ModelsA model that regresses each observation on a combination of past obser-vations with a random component is called an autoregressive model (AR)[4]. The number of historic observations considered for the current observa-tion will determine the order of the AR model. Specifically, an AR modelof order p and mean µ is given asyt = µ+ ϕ‘(yt−1 − µ) + ϕ2(yt−1 − µ) + · · ·+ ϕp(yt−p − µ) + εt. (2.6)where {εt} is an independent sequence of normally distributed random vari-ables with zero mean and constant variance σ2 and {ϕi} are constants [13].To obtain the 0-mean AR(p) model set zt = yt − µ which yieldszt = ϕ1zt−1 + ϕ2zt−2 + · · ·+ ϕpzt−p + εt, (2.7)52.1. TIME SERIES FOUNDATIONSwhere assumptions are derived from equations (2.6) and [13].Reducing the order of an AR(p) to one yields an AR(1) model. Thetypical form follows from a general AR model, where the 0-mean AR(1)model is given in [13] aszt = ϕzt−1 + εt, (2.8)and the errors are independent εt ∼ N (0, σ2ε). Under the stationary as-sumption, the initial point is assumed to be independent and follows z0 ∼N (0, σ2ε/(1− ϕ2)), while |ϕ| < 1.Following the model (2.8), the assumptions above and using linearityand independence, the expected value and variance of an AR(1) model areE[zt] = 0 and Var(zt) = σ2ε/(1−ϕ2). The derivations of these results can befound in Appendix C (see Equations (5.1) and (5.4)).The variance expression for the AR(1) model holds when |ϕ| < 1, other-wise the variance would not converge to σ2ε/(1− ϕ2) [13].It can be shown that for an AR(1) model γk = ϕkσ2z and ρk = ϕk [13].The derivations of these can be found in Appendix C (see Equations (5.5)and (5.6)).2.1.3 Simple Exponential Smoothing MethodsIn this section, exponential smoothing methods are introduced. Specifi-cally, forecasting methods based on exponential smoothing, with a focus onthose utilized in the framework proposed by the authors in [8]. This willmotivate the work in [8] and the majority of the work in this thesis.Exponential smoothing is a forecasting technique that uses weighted av-erages of past observations to produce forecasts, the weights decreasing ex-ponentially as the observations become more historical [16]. The basic formof exponential smoothing is called simple exponential smoothing, which fa-cilitates modeling and forecasting on time series data that do not show aclear trend or seasonality [16]. This class of smoothing methods has been inuse since the 1950s to help produce time series forecasts [17].The Holt-Winter Linear Trend (HWLT) method is an extension of simpleexponential smoothing, such that forecasts can be obtained when the dataexhibit seasonality and a trend [16]. The level provides the average value oftime series at time t, while the trend models the long-run movement of thelevel. The smoothing parameters for level and trend are denoted by α andβ, where 0 ≤ α, β ≤ 1. The number of steps in advance that are forecast isdenoted by h. We follow the notation and information from [11], the modelformulation is62.1. TIME SERIES FOUNDATIONSℓt = αyt + (1− α)(ℓt−1 + bt−1) (2.9a)bt = β(ℓt − ℓt−1) + (1− β)bt−1 (2.9b)yˆt+h|t = ℓt + hbt. (2.9c)Equations (2.9a) and (2.9b) represent the estimate of the level and trendat time t, respectively. Since this is an extension of simple exponentialsmoothing, (2.9a) is the weighted average of the observation yt and theone-step training forecast ℓt−1 + bt−1 [16].The initial values for the level (2.9a) and the trend (2.9b) are estimatedby minimizing the sum of squared errors (SSE) [16]. Specifically, minimizingSSE for training errors in one step by finding the values of the initial andunknown parameters [16]. The SSE from [16] isSSE =n∑i=1(yt − yˆt|t−1)2 =n∑i=1e2t . (2.10)The authors in [30] extend the HWLT method to integrate two sea-sonal components, which will allow implementation in more complex timeseries patterns. The method is called the Holt-Winter linear trend extension(HWLTE) method, which is given in [8] asyt = ℓt−1 + bt−1 + s(1)t + s(2)t + dt (2.11a)ℓt = ℓt−1 + bt−1 + αdt (2.11b)bt = bt−1 + βdt (2.11c)s(1)t = s(1)t−m1 + γ1dt, s(2)t = s(2)t−m2 + γ2dt. (2.11d)Equation (2.11a) represents two seasonal components, along with their re-spective periods denoted by m1 and m2. Equations (2.11b) and (2.11c)represent the level and trend in time t. The error dt is a random variable ofwhite noise and the set of smoothing parameters are given by {α, β, γ1, γ2}.The work in [8] integrates this method for the construction of their two newmethods to forecast complex seasonality in time series. The Holt-Winteradditive (HWA) method was introduced in [40].The HWA method was built from simple exponential smoothing andHWLT methods, so that seasonality can be captured consistently and accu-rately [16]. Specifically, the HWA method is preferred when the time series72.1. TIME SERIES FOUNDATIONSdata exhibit consistent seasonality over the given time period [16]. More-over, this refers to the seasonality remaining consistent with the level of thetime series, showing consistent patterns and amplitude over time.The HWA method and assumptions used in [16] are given byℓt = α(yt − st−m) + (1− α)(ℓt−1 + bt−1) (2.12a)bt = β(ℓt − ℓt−1) + (1− β)bt−1 (2.12b)st = γ(yt − ℓt−1 − bt−1) + (1− γ)st−m (2.12c)yˆt+h|t = ℓt + hbt + st+h−m(k+1). (2.12d)The level and trend equations for the HWA method are given by (2.12a) and(2.12b). This method incorporates a smoothing equation given by (2.12c)for the seasonal component. The seasonal component is expressed as theweighted average of the current seasonal index (yt−ℓt−1−bt−1) and the sea-sonal index m periods ago (st−m) [16]. Lastly, the level component reflectsthe weighted average of seasonally adjusted observations (yt − st−m) andthe non-seasonal forecast (ℓt−1 + bt−1), while the trend component remainsconsistent with the HWLT method [16].Similarly to the HWA method, the Holt-Winter Multiplicative (HWM)method was also introduced in [40]. This method differs from the HWAmethod as it is preferred when time series data exhibit seasonal patternsthat vary with the level of the series [16].The notation and assumptions used in [16] are given byℓt = αytst−m+ (1− α)(ℓt−1 + bt−1) (2.13a)bt = β(ℓt − ℓt−1) + (1− β)bt−1 (2.13b)st = γyt(ℓt−1 + bt−1)+ (1− γ)st−m (2.13c)yˆt+h|t = (ℓt + hbt)st+h−m(k+1). (2.13d)The level equation given by (2.13a) is the weighted average of seasonallyadjusted observation yt/st−m and non-seasonal forecast ℓt−1+bt−1. Further-more, (2.13c) represents the smoothing equation for the seasonal component,which is the weighted average of the current seasonal index yt/(ℓt−1 + bt−1)and the seasonal indexm periods ago st−m [16]. Lastly, the trend component(2.13b) remains consistent with the HWA method [16].82.1. TIME SERIES FOUNDATIONS2.1.4 Modified Models: TBATS and BATS ModelsTo extend established exponential smoothing methods, [8] propose meth-ods that can handle complex seasonal time series using exponential smooth-ing. The first model is the Box-Cox transform, ARMA errors, trend and sea-sonal components (BATS) model. This is an extension of the HWLTE givenby (2.11a)-(2.11d). The second model is the trigonometric BATS model(TBATS), which combines seasonal components with the BATS model us-ing trigonometric representations based on Fourier series from [39]. Thesetwo models from [8] aim to provide accurate and consistent forecasts forcomplex seasonal time series while addressing issues such as nonlinearity.The BATS model extends (2.11a)-(2.11d) [8], where the form is given asy(ω)t ={yωt −1ω if ω ̸= 0,log(yt) if ω = 0,(2.14a)y(ω)t = ℓt−1 + ϕbt−1 +T∑i=1s(i)t−mi + dt, (2.14b)ℓt = ℓt−1 + ϕbt−1 + αdt, (2.14c)bt = (1− ϕ)b+ ϕbt−1 + βdt, (2.14d)s(i)t = s(i)t−mi + γidt, (2.14e)dt =p∑i=1ψidt−i +q∑i=1θiεt−i + εt. (2.14f)The Box-Cox-transformed observations are given by (2.14a), where ω ̸= 0yields a power transformation, while ω = 0 employs a logarithmic transfor-mation. This limits the study to positive time series, which is quite commonin applications [8]. The damping parameter is given by ϕ, where the authorsfollow the work on damped trends from [11]. The seasonal periods are givenby {m1, . . . ,mT }, b is the long-term trend, bt is the short-term trend, dtdenotes an ARMA(p, q) process, and εt ∼ N(0, σ2). The level, trend andseasonal smoothing parameters are given by {α, β, γi} for i ∈ {1, . . . , T}.The addition of the long-term trend b allows the predictions of bt → binstead of zero [8]. Including the damping parameter ϕ in (2.14b), (2.14c),and (2.14d) ensures consistency with the work completed in [11]. The set ofparameters for a BATS model is given as BATS{ω, ϕ, {p, q},m1, . . . ,mT },which denotes the use of Box-Cox transformations, orders of ARMA errors,and seasonal periods.92.2. STATE SPACE MODELSThe TBATS model extends the BATS model given by (2.14a)-(2.14f) byaltering the seasonal components ass(i)t =ki∑j=1s(i)j,t , (2.15a)s(i)j,t = s(i)j,t−1 cos(λ(i)j ) + s∗(i)j,t−1 sin(λ(i)j ) + γ(i)1 dt, (2.15b)s∗(i)j,t = −s(i)j,t−1 sin(λ(i)j ) + s∗(i)j,t−1 cos(λ(i)j ) + γ(i)2 dt, (2.15c)as given by [8]. This flexible approach to seasonal components is based on theFourier series of [39]. The smoothing equation (2.15a) is substituted in theBATS model for (2.14e). Furthermore, (2.15b) represents the stochastic levelof the ith seasonal component, while (2.15c) denotes the stochastic growthof the ith seasonal component. The set of seasonal smoothing parametersis given by γ(i)1 and γ(i)2 for i ∈ {1, . . . , T}.The angular frequency is given by λ(i)j = 2πj/mi, where mi representsthe ith seasonal component. The angular frequency determines how theith seasonal component oscillates. The authors in [8] denote the number ofharmonics for the ith seasonal component by ki. This provides the set ofparameters for a TBATS model, which is given asTBATS{ω, ϕ, {p, q}, {m1, k1}, . . . , {mT , kT }}.2.2 State Space ModelsState space models are a common approach to time series modeling,forecasting, and clustering. The flexible framework provided by state spacemodels through the observed and state (latent) system equations can capturetime-series patterns such as seasonality and trend [26]. For time series,state space models model serial dependence between observations, enablinga framework that can account for temporal structure in time series data[3]. For example, state space models can be useful in applications includingecological processes, population dynamics, economics, and more [3].Applications where the time series formulation is known or specified en-able assumptions on state space models to be made, which can determinemodel selection criteria, likelihood computations, confidence intervals, andmore [26]. The flexible framework provided from state space models arisesfrom the observed equation for the observations produced from the state(latent) equation, that models the hidden and true underlying dynamics of102.2. STATE SPACE MODELSthe system [3]. The implementation of state space models provides the ad-vantages of a flexible approach, an interpretable framework, handles missingvalues, and model selection criteria [26, 35].2.2.1 Linear State Space ModelsA typical way to define state space models can be found in [19], whereit is assumed that yt is a dy-variate time series and xt is a dx-dimensionalunobservable vector [19]. The observed equation is denoted by yt, whilethe state (latent) equation is denoted by xt. A linear state space modelis constructed from a linear relationship between the observed and stateequations in the system; A typical formulation for a linear state space modelfor a time series yt can be given asyt = Htxt + vt (2.16a)xt = Ftxt−1 +Gtwt. (2.16b)It is assumed that wt is a dm-dimensional white noise vector with meanzero and variance-covariance matrix Qt and vt is a dy-dimensional Gaussianwhite noise vector with zero mean and variance-covariance matrix Rt. Fur-thermore, the transition matrix Ft, the observation matrices Ht, and Gt are(dx × dx), (dx × dm), and (dy × dx). The assumptions are from [19].To provide an example of linear state space models in time series, con-sider the 0-mean moving average model of order 2 (MA(2)). Following thenotation in [13], the model can be given aszt = θ1εt−1 + θ2εt−2 + εt. (2.17)Here {θ1, θ2} are constants and εt ∼ N(0, σ2ε) are independent random vari-ables [13]. Let the state vector be xt = (εt, εt−1, εt−2)⊤, then a state spacerepresentation for the 0-mean MA(2) model can be given as:xt =0 0 01 0 00 1 0︸ ︷︷ ︸Ftεt−1εt−2εt−3︸ ︷︷ ︸xt−1+100︸︷︷︸Gtεt︸︷︷︸wt, (2.18)yt =[1 θ1 θ2]︸ ︷︷ ︸Ht εtεt−1εt−2︸ ︷︷ ︸xt+ 0︸︷︷︸vt. (2.19)112.2. STATE SPACE MODELS2.2.2 Linear Gaussian State Space ModelsA powerful special case of a linear state space model is the linear Gaus-sian state space model (LGSSM). Consider the observed and state variablesgiven by Y ={yt}Tt=1, X ={xt}Tt=1. For LGSSMs, it is assumed thatboth X and Y follow multivariate Gaussian distributions. The form andassumptions from [35] are given asxt = Axt−1 + wt (2.20a)yt = Cxt + vt (2.20b)x1 = µ+ u. (2.20c)It is assumed that yt is a dy-variate time series while xt is a dx-dimensionalunobservable vector. The noise terms are defined as wt ∼ N (0,Γ), vt ∼N (0,Σ), and u ∼ N (0, P ). Γ, Σ, and P are (dx×dx), (dy×dy), and (dx×dx)variance-covariance matrices respectively. Furthermore, A is known as the(dx× dx) transition matrix, C is the (dy× dx) observed matrix, and µ is theinitial dx-dimensional vector.The set of parameters given by θ ={A,Γ, C,Σ, µ, P}is typically ofinterest for a given LGSSM. The observation matrix C is set to C = [1 . . . 1]in [35] and is not estimated since the LGSSM representation is not uniqueotherwise for a given set of observed variables. A detailed description of theestimation for θ for an LGSSM can be found in [35].To provide an example of a time series model as an LGSSM, considerthe AR(2) model with the notation of [13], where εt ∼ N(0, σ2e) andzt = ϕ1zt−1 + ϕ2zt−2 + εt. (2.21)Now set the state vector as xt = (zt, zt−1)⊤, which provides the first LGSSMrepresentation asxt =[ϕ1 ϕ21 0]︸ ︷︷ ︸A[zt−1zt−2]︸ ︷︷ ︸xt−1+[εt0], yt =[1 0]︸ ︷︷ ︸C[ztzt−1]︸ ︷︷ ︸xt+ vt. (2.22)Now change the state vector to xt = (zt − zt−1, zt−1)⊤, which gives thesecond LGSSM representation asxt =[(ϕ1 − 1) (ϕ2 + ϕ1 − 1)1 1]︸ ︷︷ ︸A[zt−1 − zt−2zt−2]︸ ︷︷ ︸xt−1+[εt0], (2.23)122.2. STATE SPACE MODELSyt =[1 1]︸ ︷︷ ︸C[zt − zt−1zt−1]︸ ︷︷ ︸xt+ vt. (2.24)The new definition of the state vector for (2.23) alters the transition matrixA, xt−1, and the observed matrix C compared to (2.22). However, the out-come of the observed variables is the same. That is, for (2.22) the observedvariables can be expressed asyt =[1 0]︸ ︷︷ ︸C[ztzt−1]︸ ︷︷ ︸xt+vt= zt + vt.For the second LGSSM representation in (2.24) the observed variables aregiven byyt =[1 1]︸ ︷︷ ︸C[zt − zt−1zt−1]︸ ︷︷ ︸xt+ vt= zt − zt−1 + zt−1 + vt= zt + vt.Both the LGSSM representations as in (2.22) and (2.23) and (2.24) assumedifferent states xt, but fail to provide unique LGSSM representations.2.2.3 State Estimation TechniquesEstimation of the state equation can introduce difficulties for a givenstate space model; Inaccurate data, noisy data, and insufficient number ofobservations make estimation of the state equation problematic at times[25]. Various estimation techniques have been constructed for state spacemodels to tackle problems in estimating the state.One of the more common techniques in linear state space models is theKalman filter and its extensions [25]. A limitation imposed by the Kalmanfilter is the linear relationship between the observed and state equations,which limits the applicability to certain systems [25]. The work in thisthesis, involving the exponential smoothing methods BATS and TBATS, isconstructed with linear state space equations, provided from [8]. Although132.2. STATE SPACE MODELSnot directly implemented, the Kalman filter is widely used in state spacemodel estimation, making it worthwhile to discuss in this work, for comple-tion and future extensions.The Kalman filter is a recursive approach to estimate the state of a statespace model. The two main steps involve prediction and updating [38]. TheKalman filter prediction step provides a prior estimate of the state at timet (xˆt|t−1), based on observations up to time t − 1, while the update stepcorrects this prediction using the new observation yt to produce the filteredestimate xˆt|t [38]. The prediction step utilizes previous knowledge on thesystem through the state equation, while the update step incorporates boththe previous estimation of the state and current observations to update stateestimates [18]. The prediction and update steps for the Kalman filter arealternatively iterated, providing state predictions and corrections to enhancethe modeling of the system dynamics over time.Implementing the Kalman filter for innovations of state space modelsfor BATS and TBATS has not been proposed, but would be viable [16].The integration of the Kalman filter for the state space innovations of theBATS and TBATS models would impose additional work to construct newsystem equations aligned with the BATS and TBATS models [16]. Detaileddescription of the Kalman filter along with equations for the iterated processfor the observed and state equations can be found in [18, 38].2.2.4 State Space Innovations: BATS and TBATSModified models such as BATS and TBATS models are linear state spacemodels that incorporate Box-Cox transformations [8]. This framework fa-cilitates the derivation of exact likelihood estimates, forecast error variance,and forecast intervals [8].The linear innovations of state space models for the BATS and TBATSmodels are provided in [8] asy(ω)t = w′xt−1 + εt (2.25a)xt = Fxt−1 + gεt. (2.25b)Equations (2.25a) and (2.25b) represent the observed and state equations,respectively. Here, w′ is the row vector, g is the column vector, F is theunobserved matrix, xt is the unobserved state vector at time t, and εt is theobserved noise.Different linear state space representations can be constructed for theBATS and TBATS models. These representations are useful, as they provide142.2. STATE SPACE MODELSa flexible framework to incorporate trend and seasonal components whilehandling issues such as missing data and detecting the underlying structureof the time series data. Examples of linear state space representations forBATS and TBATS models can be found in [8].The BATS and TBATS model selection procedure differs from the usualapproach for linear state space models, as estimates for the unknown Box-Cox parameter ω and the ARMA(p, q) process must be obtained [8]. Atypical approach to obtain one-step-ahead prediction errors and variancesis to employ the Kalman filter [8]. The Kalman filter could be viable forthe BATS and TBATS models in terms of estimation if additional equationswere added to handle non-stationary states [8]. Instead of conditioning on allstates and treating them as unknown fixed parameters, the updated Kalmanfilter can be replaced by exponential smoothing. This leads to exponentialsmoothing providing one-step-ahead prediction errors for likelihood evalua-tion [8].Estimation of unknown parameters for innovations in linear state-spacemodels for both BATS and TBATS involves calculating the conditionallikelihood of the observed data y = (y1, . . . , yn), which provides the log-likelihood. The conditional likelihood of the observed data is derived underthe assumption that εt ∼ N(0, σ2). Consequently, the density of the trans-formed time series is given as y(ω)t ∼ N(w′xt−1, σ2). The density of thetransformed data is given as:p(y(ω)|x0, ϑ, σ2) =n∏t=1p(y(ω)t |xt−1, ϑ, σ2)=n∏t=11(2πσ2)1/2exp( −12σ2(εt − 0)2)(2.26)=1(2πσ2)n/2exp( −12σ2n∑t=1ε2t). (2.27)The vector containing the Box-Cox parameter, smoothing parameters,and ARMA coefficients is denoted byϑ [8]. The Jacobian of the Box-Coxtransformation is deployed to provide the density of yt on the original scaleas:p(yt|x0, ϑ, σ2) = p(y(ω)|x0, ϑ, σ2)∣∣∣∣det(∂y(ω)t∂y)∣∣∣∣152.3. CLUSTERING TIME SERIES MODELS=1(2πσ2)n/2exp( −12σ2n∑t=1ε2t) n∏t=1yω−1t . (2.28)The resulting log-likelihood obtained from (2.28) is given in [8] asL∗(x0, ϑ) = nlog( n∑i=1ε2t)−2(ω − 1)n∑t=1log(yt), (2.29)see Equation (5.8) and Equation (5.9) in Appendix D for the full derivationdetails.2.3 Clustering Time Series ModelsClustering time series models aim to uncover underlying patterns, fea-tures, and groups of a given time series dataset. Some practical applicationsof time series clustering have led to detection in brain activities, energyconsumption patterns, and more [35]. Time series clustering typically hasthree main approaches: model-based, shape-based, and feature-based [1].Model-based approaches in time series clustering have an advantage overother methods as they can identify the correct number of clusters givenprior information [42], while also providing more accurate results comparedto other methods, provided an appropriate model is selected [35]. Further-more, a model-based approach allows for straightforward implementationsfor future predictions and simulations from the estimated parameters foreach cluster [35].Some previous work on model-based time series clustering has focusedon AR models [36], ARMA models [42], and hidden Markov models (HMM)[12]. AR and ARMA-based model approaches are limited to stationary timeseries data and struggle to capture the underlying dynamics of the evolvingtime series system [35]. To provide accurate modeling and predictions, theselected model should have a high flexibility to account for the observed andunobserved system dynamics [35]. The implementation of state space mod-els for time series clustering provides a basis for highly flexible modelling,capable of handling missing data and non-stationarity [35].2.3.1 Mixtures of Linear Gaussian State Space ModelsThe work in [35], extends LGSSMs as a mixture of linear Gaussian state-space models (MLGSSMs), allowing for a model-based clustering framework.Following [35], where Dy = {Yi}Ni=1 represents a dataset with N total time162.3. CLUSTERING TIME SERIES MODELSseries. Each time series is generated from m different LGSSMs and letDx = {Xi}Ni=1 denote the corresponding set of latent time series for Dy. LetM represent the total number of clusters, indicated by ω(1), . . . , ω(M). Letk index the LGSSM model for cluster k, which is given in [35] asxt = A(k)xt−1 + wt (2.30a)yt = C(k)xt + vt (2.30b)x1 = µ(k) + µ. (2.30c)Similar to single LGSSMs, it is assumed wt ∼ N (0,Γ(k)), vt ∼ N (0,Σ(k)),and u ∼ N (0, P (k)). Following [35], set p(k) = p(ω(k)) represents the proba-bility associated with the group k. This results in the constraint of the prob-abilities assigned to all the clusters sums to 1, equivalently∑Mk=1 p(k) = 1.Following the single LGSSMs, each cluster will have a set of LGSSMparameters. For example, θ(k) = {A(k),Γ(k), C(k),Σ(k), µ(k), P (k)} will rep-resent the parameters that the models cluster k. In turn, [35] uses thefollowing definition for the MLGSSM distributionp(Y,X|Θ) =M∑k=1p(Y,X | ω(k), θ(k))p(k). (2.31)Together, this provides the set of parameters for the entire mixture distribu-tion given by Θ = {(θ(1), p(1)), . . . , (θ(M), p(M))}. This set of parameters is ofinterest when implementing the expectation maximization (EM) algorithmfor MLGSSMs. A detailed summary of the EM algorithm for MLGSSMscan be found in [35], where the EM algorithm provides the assignment ofthe clusters for each time series along with the set of parameters for eachcluster.The MLGSSM clustering framework serves as a basis for constructingthe model-based time series clustering approach taken in this thesis onweekly wildfire occurrences in Canada. Specifically, the use of the BATSand TBATS models for each time series aligns with the modeling of eachtime series as a linear state space model [35]. Once this is complete, aclustering approach from a finite Gaussian mixture model using the EMalgorithm will be constructed on the weekly wildfire occurrences similar tothe clustering approach for MLGSSM in [35].172.4. MARKOV CHAINS2.4 Markov ChainsMarkov chains are stochastic processes that form a sequence of randomevents, where the probability of each event in the state space depends onlyon the previous state. They are widely used in physical, biological, economic,and other fields to model a variety of phenomena [9]. Formally, {Xn}n∈N iscalled a discrete-time Markov chain ifP (Xn+1 = xn+1 | Xn = xn, . . . , X0 = x0) = P (Xn+1 = xn+1 | Xn = xn)= p(i, j) (2.32)where p(i, j) denotes the transition matrix [9]. Equation (2.32) highlightsthe memoryless property of Markov chains, where the probability of the nextstate Xn+1 depends solely on the current state Xn. Although memory isrestricted to one previous state in this case, it can be extended to include afinite number, allowing for more flexible modeling.A Markov chain of finite order is denoted similarly to equation (2.32),but in this case, the future state depends on a finite number of lagged states[5]. For example, consider finite-order Markov chains (FOMCs) up to anorder of k:P (Xt+1 = xt+1 | X1, . . . , Xt) = P (Xt+1 = xt+1 | Xt) (Order 1)P (Xt+1 = xt+1 | X1, . . . , Xt) = P (Xt+1 = xt+1 | Xt−1, Xt) (Order 2)...P (Xt+1 = xt+1 | X1, . . . , Xt) = P (Xt+1 = xt+1 | Xt−k+1, . . . , Xt).(Order k)The Markov chains of varying finite orders above all have a finite memoryprocess that is constant for every state in the given state space X [24]. AMarkov chain of order k is called full if the state space is the entire X k,meaning that all future states depend on all possible length-k histories [24].FOMCs are useful in a variety of applications and fields, but have lim-ited structural richness and dimensionality issues [24]. Specifically, FOMCsoften struggle to capture meaningful and complex patterns in data with rea-sonable model selection, which grows exponentially in complexity as orderk increases [24].182.4. MARKOV CHAINS2.4.1 Variable Length Markov ChainsAn alternative to FOMCs is the variable length Markov chain (VLMC),a high-order Markov chain with a state space X from a finite categoricalset. A VLMC adaptively determines the memory length required for thetransitions [24]. In other words, it considers a subset of lagged historicalstates when transitioning to a future state, the history length varying basedon the current state in X , thus reducing the computational load of high-orderFOMCs.To introduce VLMCs formally, an example inspired by the CNFDBdataset and [24] will illustrate the typical process. Assume Xt is a sequenceof categorical outcomes drawn from a subset of the CNFDB dataset, whereX = {H,N,U} denotes the possible causes of wildfires, and the processassumes an order of 2. The structure of the transition probabilities can bedefined as:P [Xt = xt | Xt−1 = xt−1, Xt−2 = xt−2] =P [Xt = xt | Xt−1 = H], if xt−1 = HP [Xt = xt | Xt−1 = N ], if xt−1 = NP [Xt = xt | Xt−1 = U,Xt−2 ∈ {H,N}], if xt−1 = N, xt−2 ∈ {H,N}P [Xt = xt | Xt−1 = U,Xt−2 = U ], if xt−1 = U, xt−2 = U(2.33)This formulation of Xt illustrates that the transition probabilities dependon a varying number of past states, determined by the history of laggedvalues. For example, when the state Xt−1 = U , the model incorporates theprevious state at Xt−2 to determine the transition probability. In contrast,when Xt−1 = H, no additional information is needed beyond the prior state.Given the structure of Xt, the transition probabilities can be definedand denoted by πi, which is a card(X )×1 vector composed of the transitionprobabilities for each state combination of X [24]. Specifically, based onthe previous example of a second-order VLMC with X = {H,N,U}, thegiven structure, and the work in [24], these transition probabilities can beexpressed as:xt−1t−2 H∗ N∗ U [HN ] UUP [Xt | · ] π1 π2 π3 π4Following the notation in [24], for example, H∗ denotes Xt−1 = H and192.4. MARKOV CHAINSXt−2 = arbitrary, while U [HN ] denotes Xt−1 = U and Xt−2 ∈ {H,N}.The probability distributions associated with the state transitions from Xare indicated by π1, π2, π3, π4, where the components of πj sum to one [24].For example, let π3 = (1/3, 1/3, 1/3), which would define the following:P [Xt = H | Xt−1 = U,Xt−2 ∈ {H,N}] = 1/3P [Xt = N | Xt−1 = U,Xt−2 ∈ {H,N}] = 1/3P [Xt = U | Xt−1 = U,Xt−2 ∈ {H,N}] = 1/3.The previous example has four distinct transition probability vectors πjin comparison to a generic full-order 2 FOMC, which would have card(X )2 =9. This illustrates how constructing a VLMC can reduce the number oftransition probability vectors required by varying the history of lagged valuesand using a more flexible approach to modeling the dependencies betweenstates. This results in a more compact and efficient representation comparedto traditional Markov chain models [24].2.4.2 Variable Length Memory, Context Functions, andContext TreesThe previous section served as an introduction to VLMCs through ex-amples from the CNFDB dataset. To formally define VLMCs, the variablelength memory, context functions, and context trees must be discussed.Assume (Xt)t∈Z is a stochastic process with a categorical space X , thendefinecpre : X∞ →∞⋃j=0X j ∪ X∞ (X 0 = ∅), (2.34)as the preliminary function that maps the infinite past to a possible shorterstring, representing the relevant history [24]. Moreover, the authors in [24]also denotecpre : x0−∞ 7→ x0−ℓ+1, where ℓ is defined by (2.35)ℓ = ℓ(x0−∞)= min{k : P [X1 = x1 | X0−∞ = x0−∞] (2.36)= P [X1 = x1 | X−k+10 = x−k+10 ] ∀ x1 ∈ X}, (2.37)and ℓ ≡ 0 corresponds to independence, (2.38)202.4. MARKOV CHAINSwhere cpre(·) determines which vectors are relevant from the infinite pastand is called the preliminary context function. The parameter ℓ determinesthe length of the relevant history for cpre(·) by considering the infinite pastup to the current state (x0−∞). Specifically, ℓ determines the smallest k suchthat the probability distribution for the next state, which depends on theprevious k states, is approximately equal to the distribution based on theentire history [24]. In other words, ℓ represents the variable memory lengthof the given context, which captures the minimal amount of informationfrom the past that is needed to predict the next state for the given context.Let t ∈ Z, then cpre(xt−1−∞) is termed the context process at time t [24].Additionally, let k ∈ (0,∞) be the smallest integer such thatcard(cpre(x0−∞)) = ℓ(x0−∞) ≤ k, ∀ x0−∞ ∈ X , (2.39)then k is termed the order of the VLMC. If k < ∞, the process (Xt)t∈Zis called a stationary VLMC of order k [24]. Furthermore, if the VLMC isof order k and the full projection x0−∞ 7→ x0−k+1 holds for all x0−∞, then(Xt)t∈Z is considered a full VLMC of order k [24].Given the definitions of cpre(·), we can form the standard context func-tion c(x0−∞). To illustrate this, consider the VLMC from the CNFDBdataset with X = {H,N,U} of order 2, whose structure is given in equa-tion (2.33). Using these definitions and previous results, we determine thecontext function as follows:c(x0−∞) =H, if x0 = H,x−1−∞ arbitraryN, if x0 = N, x−1−∞ arbitraryU [HN ], if x0 = U, x−1 ∈ {H,N}, x−2−∞ arbitraryUU, if x0 = U, x−1 = U, x−2−∞ arbitrary,(2.40)which is composed of 4 contexts {H,N,U [HN ], UU}, all of which are oforder 2 or lower. In summary, equation (2.40) defines the context functionc(x0−∞), which classifies the process’s history based on recent states. Thecontext function assigns a specific context to each combination of past states,capturing the relevant history for predicting the next state, and identifyingdistinct patterns in the process. Additionally, it provides structure to theVLMCs, which will be important when constructing these models [24].The graphical representation of a VLMC is given by a context tree [24].Specifically, a context tree denotes the minimal state space of the underlyingVLMC. Assume c(·) is the context function for a VLMC; then the contexttree τ is defined in [24] as212.5. WILDFIRE OCCURRENCES IN CANADAτ = τc = {w ; w = c(x0−∞), x0−∞ ∈ X∞}. (2.41)Each path w ∈ τc represents a context for the VLMC, where {w1, w2, . . . }forms the minimal state space. This structure allows c(·) to be constructedfrom τc, and vice versa.Consider the VLMC from the CNFDB dataset with X = {H,N,U} oforder 2, whose structure is given in equation (2.33) and the context functionby equation (2.40). The context tree τc is summarized below:H N UU HNFigure 2.1: Context tree τc for X = {H,N,U} from a VLMC fit fromAlberta, 2023. Each path from the root to a leaf denotes a context w.Context functions and trees are essential to define and understand thestructure of a VLMC. The context function determines which parts of thepast are relevant for prediction, while the context tree offers a clear graphicalrepresentation of this structure and the minimal state space. Together, theydefine the structure of a VLMC by the minimal state space, capturing thekey dependencies in the process without an overly complex framework.2.5 Wildfire Occurrences in CanadaWildfire occurrences are projected to increase in coming years, withprovincial variations [41]. Some factors influencing this increase in fire oc-currences are human activity and an increase in lightning-caused fires [14].In addition, unknown causes of fires typically account for the remainingcauses [33]. This places most of the causes of wildfires in Canada as human,natural, and unknown. Human and natural causes of wildfires in Canadatend to account for most fires and are evenly split [22].222.5. WILDFIRE OCCURRENCES IN CANADAWildfires in Canada can vary in terms of the negative impacts on soci-ety and the environment [33]. In addition, not all wildfires can cause majordestruction and damage to the environment and property. For example,in 2017 $137.3 million was claimed through insurance companies in BritishColumbia alone [33]. This illustrates that overtime, wildfire destruction cancompound and adhere to massive insurance claims and property damage.Furthermore, visible short-term destruction is not all that should be con-sidered, aspects such as depreciating property values, community mentalhealth issues, and business revenue should be included in total costs [10].In addition, the most destructive wildfires in Canadian history tend to beattributed to human causes, where the major structural damage, humandeaths, and injuries increase [22].Wildfire management in Canada can be broken down into the followingmanagement groups: prevention, mitigation, response, preparedness, andrecovery [23]. To combat the spread and destruction of wildfires in Canada,these management groups should be implemented. Specifically, to optimizeresource allocation, forecast accuracy and duration of fire occurrences areessential [41].This presents a problem for studying and improving current forecastmethods for wildfire occurrences in British Columbia and throughout Canada.Providing reliable and consistent wildfire occurrence forecasts can be vitalto optimizing wildfire risk groups. Specifically, forecasts can provide insightinto the number of wildfire occurrences across Canada for a given period oftime or location. In addition, studying the causes of wildfires in Canadashould also be of concern. Specifically, the ability to forecast and simulatethe causes of wildfires can help identify trends in all regions of Canada. Inturn, this may help provide insight into future causes and help implementpreventative measures across Canada.A current method in place in British Columbia provides reliable 7-dayforecasts for weekly fire counts, providing a method to identify future pat-terns. The framework to construct will be concerned with determining whenfluctuations might occur, rather than how many wildfires will occur in aprovince over the given period. In addition, the new method will be used toidentify the time series class of weekly wildfire occurrences for an upcomingseason.2.5.1 Wildfire Occurrences DatasetThe CNFDB dataset records forest fire information from provincial andterritorial fire management agencies throughout Canada [6]. To account for232.5. WILDFIRE OCCURRENCES IN CANADAthe variation in documentation from the CNFDB dataset by region, forestfires are documented as a subset of wildfires in this thesis, which is true ingeneral. The CNFDB dataset records 436,564 fires in all Canadian provincesfrom 1930 to 2023. Some of the attributes included are fire ID, location, startand end dates, size and cause [6]. For this thesis, a subset of variables willbe used to model and forecast wildfire occurrences and causes. A notablelist of the variables from the CNFDB dataset can be found in Appendix ATable 5.1.The causation of every wildfire is classified into five categories: human-caused (H), natural-caused (N), unknown cause (U), human-prescribed burn(HPB), and reburn (RE). Of the 436,564 wildfires, the distribution is asfollows: 235,898 human-caused (H), 189,625 natural-caused (N), 10,626 un-known (U), 341 human-prescribed burns (HPB), and 74 reburns (RE) [6]. Toaddress class imbalances, this analysis is restricted to causes X = {H,N,U},which account for the majority of wildfires and tend to be more destructive[22].The typical time series pattern for wildfire occurrences in Canada over acalendar year illustrates two large influxes of wildfire occurrences. The firstincrease occurs around early spring and the second around middle to latesummer. These two increases in wildfires are also associated with some ofthe most destructive wildfires, where characteristics such as wind, droughtperiods, and more drive this increase [22].To illustrate this seasonal pattern, Figure 2.2 shows the aggregatedweekly wildfire occurrences in British Columbia in 2003. This figure providesa generalization of what to expect, where the first increase in the number ofwildfires occurs around late May, and the largest increase appears aroundearly August. In addition, the beginning and end of each calendar yearshould also be a reduced number of wildfires compared to other months.242.5. WILDFIRE OCCURRENCES IN CANADAFigure 2.2: Weekly wildfire occurrences in British Columbia in 2003, withpeaks in late May and early July, with minimal activity during wintermonths. Data sourced from Natural Resources Canada CNFDB dataset.Extending weekly wildfire counts by province and territory over multipleyears yields a different result. Specifically, multiple individual time seriesare conjoined into an extended version that exhibits patterns of non-integerfrequency, complex seasonality, and more. This new time series of weeklywildfire counts will now be over multiple calendar years, where the sametypical pattern as in Figure 2.2 is observed over time.To visualize this extension of wildfire counts over multiple calendar years,consider Figure 2.3. This shows the aggregated number of weekly wildfiresin British Columbia from 1950 to 2023. Each year exhibits similar patternsto Figure 2.2, where the long-term seasonality and frequency differ. Eachregion within Canada should exhibit similar long-term patterns, where thecalender years within these ranges might differ.252.5. WILDFIRE OCCURRENCES IN CANADAFigure 2.3: Weekly wildfire occurrences in British Columbia from 1950 to2023. Data sourced from Natural Resources Canada CNFDB dataset.26Chapter 3Methodology3.1 Wildfire Occurrence ModelingThe purpose of this section is to describe the process used to prepare andmodel data on the occurrence of wildfires using the TBATS and BATS mod-els. Specifically, weekly counts for all provinces and territories were obtainedfrom the CNFDB dataset [6]. This procedure aims to establish a methodfor analyzing regional weekly counts to produce forecasts, simulations, and aclustering method that forms a basis for performing model-based clusteringto identify an upcoming wildfire season.The outline prepares the data for exponential smoothing models. Weeklywildfire counts were aggregated for each province and territory. These timeseries were either combined across years for forecasting or kept separate byyear for clustering. Each region was independently analyzed by specifyingits identifier, as uniformity between provinces is not assumed. Within eachprovince, the number of wildfires was counted per week for a given calendaryear. Initially, each year is treated as its own time series, but sequential yearswill eventually be joined to support intermediate- to long-term forecasting.In addition, a similar framework is applied for the causes of wildfires, wherethe investigation focuses on X = {H,N,U} and the performance of VLMCson various categorical time series. Specifically, each wildfire has a docu-mented cause, instead of aggregating weekly wildfire counts for a regionover a specific interval.The weekly wildfire occurrence time series is denoted yt, and the modelspecifications are discussed when introducing the BATS and TBATS models.Each year was treated as an individual time series and aggregated overmultiple years or left separately.3.1.1 Dataset Description and ChallengesThe CNFDB dataset was preprocessed before implementation and ap-plication. Specifically, unnecessary variables such as fire identification num-ber, national parks, and more were removed. Most of the variables removed273.1. WILDFIRE OCCURRENCE MODELINGwere not relevant to this work. Before aggregating wildfire counts and oc-currences throughout Canada, dates were converted into day, month, andyear for analysis. Once the CNFDB data set was cleaned and formattedproperly, aggregation could be deployed on the number of occurrences andcauses of wildfires by province or territory.3.1.2 Model Selection and FittingTo preprocess the time series, the inputs were transformed to log(1+ yt)to handle zeros. The resulting series were fitted using the functions bats()and tbats() from the forecast package in R by [15], as described in [8].These functions accept fixed or automatic-model specifications, includingtrend, damping, seasonal periods, harmonics, and ARMA errors.The selection of the model was based on the AIC criterion, which isautomated by the functions bats and tbats() functions from the forecastpackage in R [15]. The AIC implemented for the BATS and TBATS modelsisAIC = L∗(xˆ0, ϑˆ) + 2K, (3.1)where L∗(xˆ0, ϑˆ) denotes the log-likelihood of (2.29) for the respective model,and K is the number of parameters in ϑ [8]. The BATS or TBATS modelwill have the lowest AIC value for the time series data inputted.To complete the model selection process, it is necessary to determinethe number of harmonics ki for the ith seasonal component of the TBATSmodels and the ARMA orders p and q for both models. To determine thenumber of harmonics for a TBATS model, begin by de-trending the firstfew seasons of the data and then approximating the de-trended data using∑Ti=1∑kij=1(a(i)j cos(λ(i)j t) + b(i)j sin(λ(i)j t)) [8]. Start with a single harmonicand gradually add more using F-tests for each addition, where k∗i denotesthe number of significant harmonics with p < 0.001 [8]. The model is then fitwith ki = k∗i , and the corresponding AIC is computed. One seasonal compo-nent is considered at a time by re-fitting to the training set while graduallyincreasing ki, holding all constant until the minimum AIC is achieved.The first step in determining the ARMA orders p and q in [8] requirescompleting the process described above, assuming that there are no ARMAerror terms. Once this is complete, a model without ARMA errors isachieved, and the automatic ARIMA algorithm is applied to the residu-als of the model. This approximates the ARMA orders p and q, as indicatedin [8]. The chosen model is then re-fit to the data with the ARMA(p, q)283.1. WILDFIRE OCCURRENCE MODELINGcomponent. If the resulting AIC value for the model that includes ARMAerrors is lower than that of the model that does not have ARMA errors,then the model with ARMA errors is selected.Each fitted model produces a list object that supports the forecastingand simulation of the forecast package in R [15]. The following algorithmsummarizes the data pre-processing and automatic BATS and TBATS fit-ting procedure. Each model fit contains multiple components; details of theBATS and TBATS output structures can be found in Table 5.2.Algorithm 1 Preparing and Fitting TBATS and BATS Models to WeeklyTime Series DataRequire: Data frame D that includes a date column d and region columnr; model type ("BATS" or "TBATS")Ensure: List of fitted model objects {model1, . . . , modeln}1: Identify all distinct region and year pairs {(r1, y1), . . . , (rn, yn)} in D2: For every pair (ri, yi), i = 1, . . . , n:3: Subset D where r = ri and year = yi4: Convert column d to Date5: Aggregate wildfire counts by week to obtain Yi6: Convert Yi to time series object7: Apply transformation: Zi ← log(1 + Yi)8: If model type = "BATS":modeli ← bats(Zi, seasonal.periods = 52)9: Else if model type = "TBATS":modeli ← tbats(Zi, seasonal.periods = 52.18)10: return {model1, . . . , modeln}3.1.3 Residual Diagnostics: TBATSThe majority of the reported work in this thesis that involves modeling,forecasting, and clustering focuses on TBATS models rather than BATSmodels in BC. Residual diagnostics for BC were explored for the years 1950-2023 of weekly wildfire occurrences for the TBATS models, which will beexamined further in the clustering results.Similarly to Algorithm (1), the CNFDB dataset was subset for BC where74 time series each with 52 observations were examined for weekly wildfirecounts. The residual diagnostics were completed to examine how well theTBATS models could capture the temporal structure of each time series.The residual plots were inspected for each year for white noise with no clear293.1. WILDFIRE OCCURRENCE MODELINGtrend or pattern. Autocorrelation function (ACF) plots were produced andexamined to see if the residuals still had serial dependence upon fitting.Finally, the Box-ljung test was also applied to each calendar year to formallytest for autocorrelation in the residuals.Across the residual diagnostics for BC from 1950-2023, the majorityof the years had p-values larger than 0.05 in the box test, indicating thatthere was no significant autocorrelation. Moreover, some years producedp-values less than 0.05 indicating some autocorrelation remaining in theresiduals. Furthermore, most of the residual and ACF plots illustrated minorto no significant autocorrelation for each year. Although not every yearprovided optimal residual diagnostics, most demonstrated minimal to nosignificant autocorrelation, indicating that TBATS models are appropriatefor modeling, forecasting, and clustering time series data in BC from 1950-2023.3.1.4 Forecasting and Simulating Wildfire Time SeriesOnce Algorithm (1) is complete, the resulting model list can be used toforecast and simulate weekly wildfire occurrences in Canada. The imple-mentation of forecasting and simulation is automated with forecast() andsimulate() from the forecast package in R [8, 15]. Forecasting and simu-lation can be performed in terms of the desired number of weeks forecastedor simulated, whether a single week or multiple years are desired; both arestraightforward implementations due to the forecast package [8, 15].To assess performance and generalization, time series data from weeklywildfire counts were split into training and testing sets for each region. Thisenabled out-of-sample evaluation via forecasts for both BATS and TBATSmodels. Simulated wildfire counts were also generated to study uncertaintyin future behavior. The forecast results were assessed for each region sep-arately, where a comparison could be made with the true counts with theestimated counts over all weeks considered. Time series plots for both BATSand TBATS can be obtained to visualize the forecasting and simulation re-sults for a given region. Simulations are also used to assess the consistencyand stabilization of estimates over a longer period to validate the robustnessand stability of the selected models. Finally, metrics such as mean squarederror (MSE) and mean absolute error (MAE) were compared for BATS andTBATS model fits to study the effectiveness of both models in terms ofout-of-sample accuracy.303.2. CLUSTERING WILDFIRE TIME SERIES3.2 Clustering Wildfire Time SeriesTransitioning to the clustering framework, each time series representsweekly wildfire counts for a calendar year from a specific region. The goalis to cluster similar years based on extracted features from the given expo-nential smoothing model fit. Model-based clustering is selected so that eachcluster has a representative time series model for weekly wildfire counts ina region, where prediction and inference can be made from the fitted modelfor each cluster [35].Data preprocessing is the same as in the forecasting procedure, whereAlgorithm (1) can summarize the initial steps. Specifically, when shifting tothe clustering framework, each time series object passed into Algorithm (1)contains the weekly wildfire counts for the given region for a calendar year.This results in a list of fits of the TBATS model from weekly wildfire countsover all the recorded years for each province or territory.The model-based clustering approach is inspired by [35], where similari-ties are the usage of linear state space models for each time series, while thedifference is the specific models used and the clustering execution. Specifi-cally, the work in [35] focuses on LGSSMs for each time series, while in thisthesis each weekly wildfire occurrence time series is modeled as a BATS orTBATS model. Furthermore, instead of assuming a finite mixture of BATSor TBATS models similar to the MLGSSM in [35], we extract the featuresof each BATS and TBATS fit and perform model-based clustering from aGaussian finite mixture model. The framework in this thesis and the workin [35] rely on the EM algorithm to provide cluster assignments and param-eters for each cluster, a detailed description of the EM algorithm in thisthesis is given in the subsequent sections.3.2.1 Clustering Method and Feature ExtractionOnce each time series for weekly wildfire counts over a calendar year hasbeen fit to a TBATS model, the clustering procedure can implemented. Thedata is preprocessed with Algorithm (1) and the features of each TBATS fitare extracted. To account for varying scales for all the appropriate featuresand the clustering algorithm selected, each entry is scaled to have zero meanand unit variance. Once this is complete, model-based clustering can bedeployed on all available time series. Each time series of weekly wildfireoccurrences over a calendar year is then assigned to a cluster through the EMalgorithm, which provides the cluster assignment along with the parametersof the Gaussian finite mixture model for each cluster. This can be automated313.2. CLUSTERING WILDFIRE TIME SERIESusing Mclust() in R on the scaled feature matrix [28].Following [31], a representative time series for each cluster is found byaveraging the weekly fire counts for the cluster. Once this has been com-pleted, fit a TBATS model to the representative time series model of weeklywildfire occurrences. This allows the model-based clustering to have a rep-resentative model that can be used to simulate or forecast the occurrenceof wildfires across various regions throughout Canada over a calendar year.The following Algorithm (2) describes the steps taken in detail, while Fig-ure 3.1 provides a broad visual summary of the constructed framework.Algorithm 2 Model-Based Clustering of Time Series Using Fitted TBATSModelsRequire: List of fitted TBATS model objects {model1, . . . , modeln} fromAlgorithm (1)Ensure: Cluster assignments {c1, . . . , cn} and representative models{modelrep1 , . . . , modelrepK }1: Extract scalar features from each modeli:AIC, residuals, fitted values, variancesmoothing parameters {α, γ1, γ2}number of harmonics {ki}AR and MA parameters {p, q}2: Form feature matrix F ∈ Rn×m3: Standardize F column-wise with zero mean and unit variance4: Apply model-based clustering: Mclust(F )5: Obtain cluster assignments {c1, . . . , cn}6: For each cluster j = 1, . . . ,K:7: Identify indices: Ij = {i : ci = j}8: Extract all original series {Yi}i∈Ij9: Average weekly counts for all series for {Yi}i∈Ij10: Obtain representative series: Y¯j ← rowMeans({Yi})11: Fit TBATS model: modelrepj ← tbats(Y¯j , seasonal.periods = 52.18)12: return {ci}ni=1, {modelrepj }Kj=1323.2. CLUSTERING WILDFIRE TIME SERIES ModelModelModelModelClusteringTime Series DatasetTBATSTBATSTBATSTBATSFigure 3.1: Model-based clustering illustration with varying time seriespatterns modeled by TBATS models on each cluster representation.3.2.2 Mclust and the EM AlgorithmThe time series clustering procedure described in Algorithm (2), is amodel-based approach using the package mclust in R [28]. Specifically, thefunction Mclust() is used to perform model-based clustering on a featurematrix for weekly wildfire occurrences in a region over multiple years, provid-ing classification and density estimation results. The model-based approachoffered from Mclust(), allows for classification and density estimation re-sults from the inputted weekly wildfire occurrence data [27].Following Algorithm (2), the sample of n observations z = {z1, . . . , zn}denotes each calender year of weekly wildfire counts, while zi = {zi1, . . . , zim}denotes the set of features clustered for each calender year. Following [27],the Gaussian finite mixture model for the K clusters is given by333.2. CLUSTERING WILDFIRE TIME SERIESf(zi, φ) =K∑j=1πjfj(zi; θj), (3.2)the jth component density of zi is fi(zi; θj) for the parameter vector θjon the ith observation zi [27]. The mixing proportions are {π1, . . . , πK}(πj > 0,∑Kj=1 πj = 1) and the model parameters {θ1, . . . , θK} compose theset of parameters for the given mixture model, denoted byφ = {π1, . . . , πK , θ1, . . . , θK} [27]. The goal of Mclust() is to estimate φfrom the observation set z = {z1, . . . , zn}. This in turn will provide mixingproportions and a set of parameters for each zi, which identifies the mostlikely cluster assignment and a set of parameters θj . Since Mclust() isfrom a Gaussian finite mixture model, it is assumed that every componentfj(zi; θj) ∼ N(µj ,Σj) [27]. The mean vector µj determines the center of eachcluster and Σj determines the volume, shape, and orientation; a detaileddescription of the varying covariance structures of Σj can be found in [27].The model parameters φ are unknown for the fixed K clusters, whichmust be estimated from the input data [27]. Typically, the log-likelihood forEquation (3.2) would be derived and denoted byℓ(φ; z1, . . . , zn) =n∑i=1log(f(xi;φ)). (3.3)Maximization of the logarithmic likelihood in Equation (3.3) is not straight-forward [27]. To overcome this, Mclust() uses the EM algorithm to obtainthe MLE estimator for the Gaussian finite mixture model [27, 28]. Typi-cally, the EM algorithm iterates alternatively between the expectation step(E-step) and the maximization step (M-step) before convergence. A typicaloverview of the EM algorithm with the framework for this work is foundbelow.1. Initialization: Set initial values for φ(0) = {π(0)1 , . . . , π(0)K , θ(0)1 , . . . , θ(0)K },typically via hierarchical clustering [27].2. E-step: Compute posterior probabilities for each observation zi be-longing to cluster j using the current parameters φ(t) [37]P (cj | zi;φ(t)) =π(t)j fj(zi; θ(t)j )∑Kk=1 π(t)k fk(zi; θ(t)k ),where i = 1, . . . , n indexes observations and k = 1, . . . ,K index clus-ters.343.2. CLUSTERING WILDFIRE TIME SERIES3. M-step: Update the parameters using the current posterior probabil-ities as follows from [37]π(t+1)j =1nn∑i=1P (cj | zi;φ(t)), (3.4)µ(t+1)j =∑ni=1 P (cj | zi;φ(t))zi∑ni=1 P (cj | zi;φ(t)), (3.5)Σ(t+1)j =∑ni=1 P (cj | zi;φ(t))(zi − µ(t+1)j )(zi − µ(t+1)j )⊤∑ni=1 P (cj | zi;φ(t)). (3.6)4. Convergence: Compute the log-likelihoodℓ(φ(t+1); z1, . . . , zn) =n∑i=1log(K∑j=1π(t+1)j fj(zi; θ(t+1)j ))If the change from the previous iteration is below a threshold, termi-nate the algorithm, otherwise return to step 2 for iteration t+ 1.The final output provides the estimated parameters φˆ, posterior probabilities(cj | zi), and cluster assignments {c1, . . . , cn} for each observation zi, whichis given in the output of Algorithm (2).3.2.3 Fitting VLMCsFitting a VLMC involves implementing the tree-structured context algo-rithm, where variable-length memory is represented by an estimated contexttree [24]. Specifically, tree-structured models are fitted with each terminalnode representing x ∈ X , along with the associated transition probabilities.As noted in [24], the context algorithm grows a large tree and prunes itusing a cutoff parameter K. This parameter is used to compare a tree withits subtree by cutting a terminal node and evaluating the difference in de-viation between the two trees [24]. The value of K determines whether theVLMC fit results in a small context tree (large K) or a large context tree(small K). To determine the cutoff value for a VLMC fit, [24] recommendsdefining the cutoff asK = K(α) =12χ2ν,α, ν = card(X )− 1, (3.7)353.2. CLUSTERING WILDFIRE TIME SERIESwhere α is the specified significance level. The authors in [24] suggest usingthe K values derived from (3.7), depending on the objectives of the study.This process can be automated with a prespecified cutoff parameter K in Rusing the vlmc() function from the VLMC package [20].Once a VLMC is fitted, the model output includes the maximum memorylength, the number of unique contexts, and transition probability vectorsPˆ (xt | cˆ(xt−1, xt−2, . . . )) for each derived context. These outputs summarizethe underlying probabilistic structure and can be used to assess the fit ofthe model through the AIC or to forecast and simulate categorical statesof X based on observed context histories. This will be extended to clusterassignments from the model-based clustering approach in Algorithm (2) toidentify the time series class of the number of wildfire occurrences for anupcoming time period.36Chapter 4Applications4.1 Simulation StudiesFor the first part of the application chapter, the focus will be on theoccurrence of wildfires. Specifically, the first part will focus on a simulationstudy on TBATS models, where the consistency and reliability of simulatingfrom a single or multiple TBATS model fits will be explored. In addition,long-term simulations of weekly wildfire occurrences are examined in termsof simulating the overall pattern exhibited from the weekly wildfire occur-rences observed over multiple years.The second part of the simulation studies will look at the weekly wildfireoccurrence data for the clustering procedure. Specifically, an illustrativeexample of VLMCs on simulated clustering labels is presented in preparationfor the clustering results in subsequent sections. All applications in thissection are produced from the CNFDB dataset [6].4.1.1 TBATS Simulations for Wildfire DataIn this application, the focus is on the TBATS models and the modelingproperties for weekly wildfire counts over a calendar year and multiple years,serving as a basis for the subsequent clustering results.Specifically, Figure 2.2 is extended to include a TBATS simulation todemonstrate simulating a TBATS fit over a single calendar year, this can befound in Figure 4.1. The simulated weekly wildfire occurrences in BritishColumbia in 2003 align with the observed values. The increase or decreasein the number of fires over a given time period is captured rather than theamount. This application motivates the extension of increasing the numberof simulations and the number of years considered for weekly wildfire countsin Canada.374.1. SIMULATION STUDIESFigure 4.1: Weekly fire occurrences in BC in 2003. A TBATS model wasfit and then used to simulate wildfire occurrences for comparison.Building upon the previous single simulation from a TBATS model fit towildfire data in British Columbia in 2003, consider the full range of weeklywildfire counts in the province from 1998 to 2024. Figure 5.13 summarizesthe weekly wildfire counts observed during this period.To study simulations based in British Columbia between 1998-2023, splitthe data into a training set (1998-2017) for model selection and a testingset (2018-2023). TBATS models are fitted using a rolling window approach.For example, the first TBATS model will be trained in the years 1998-2017to simulate weekly wildfire counts for 2018. The training set will then beupdated to include the years 1998-2018 to simulate from a TBATS model for2019. This rolling window approach is completed once all years have beenindividually simulated from a TBATS fit. This sequential approach avoidssimulating the entire test period at once for a more realistic forecastingapproach.384.1. SIMULATION STUDIESFigure 4.2: Weekly wildfire counts in BC from 2018–2023 (black), overlaidwith the TBATS simulation mean (red) for the test period. Based on onesimulation trial.To extend this simulation study for TBATS models in British Columbia,the number of simulations is increased to examine how the simulated weeklywildfire occurrences stabilize, and to assess whether consistent and reliablesimulations can be obtained. The previous rolling window simulation ap-proach is conducted for the mean number of simulations. A 95% confidenceinterval will also be computed for the number of simulations used and added.The first simulation trial uses 5 iterations (Figure 4.3), while the sec-ond increases this to 50 (Figure 4.4) based on TBATS fits in the trainingset (1998-2017). For both trials, the average number of simulated weeklywildfire counts is calculated and can be compared with the observed valuesthroughout the test period. Increasing the number of simulations from 5to 50 smooths the simulated weekly wildfire counts during the test period(2018–2023). In all cases, the typical seasonal pattern is retained, where thetwo major increases typically occurring in late May and July are observed.As the number of simulations increases, smaller fluctuations are smoothedout. Individual results for the test years (2018-2023) for both simulationtrials can be found in Figures 5.1-5.6 and Figures 5.7-5.12.394.1. SIMULATION STUDIESFigure 4.3: Weekly wildfire counts in BC from 2018–2023 (black), overlaidwith the TBATS simulation mean (red) for the test period. Based on 5simulation trials with 95% confidence interval included.Figure 4.4: Weekly wildfire counts in BC from 2018–2023 (black), overlaidwith the TBATS simulation mean (red) for the test period. Based on 50simulation trials with 95% confidence interval included.404.1. SIMULATION STUDIESTo conclude this simulation study, the case with 500 simulation iterationsis also used to construct a 95% confidence interval for the simulated weeklywildfire counts. This is shown in Figure 4.5.Figure 4.5: Weekly wildfire counts in BC from 2018–2023 (black), overlaidwith the TBATS simulation mean (red) for the test period. Based on 500simulation trials with 95% confidence interval included.4.1.2 VLMC Wildfire Occurrence DataBefore entering the model-based clustering results for weekly wildfireoccurrences, a simulated example is presented to model categorical timeseries with VLMCs. Specifically, from Algorithm (2), the clustering resultswill provide a finite set of cluster assignments for each calendar year of theoccurrence of wildfires considered. Ultimately, this output can be extendedto Markov chain models and extensions to provide information on the mostlikely class or cluster of weekly wildfire occurrences over the next calendaryear. This application aims to introduce VLMCs while also motivating themodel-based clustering approach discussed in subsequent sections.For illustrative purposes, the simulated clustering labels will be restrictedto five outcomes to produce an interpretable and realistic context tree and414.1. SIMULATION STUDIESfunction. Following [24], X = {1, 2, 3, 4, 5} = {c1, c2, c3, c4, c5} will denotethe five cluster labels. There will be 500 observations, all having a cluster as-signment from X . The random distribution of outcomes for this applicationis X = {c1, c2, c3, c4, c5} = {94, 97, 106, 98, 105}.To determine K, assume that α = 0.05. Given that ν = card(X )−1 = 4,we calculate K = 12χ24,0.05 ≈ 4.74, and use K = 4 as the cutoff to producean example that is not too simple, resulting in a small context tree and astraightforward context function. Fitting this subset of simulated clusterassignments as a VLMC with the specified criteria results in the followingmodel summary:Table 4.1: VLMC model summary of 500 simulated cluster assignments.Attribute ValueCategorical Responses {c1, c2, c3, c4, c5}Observations (n) 500Cutoff Value (K) 4.5Order 3Contexts 8The previous VLMC fit of the 500 simulated cluster assignments forweekly wildfire occurrences can also provide the context function and tree.The output and corresponding context function and tree are given below.[x]-( 103 96 91 105 104| 499)+--[2]-( 21 15 23 23 14| 96)| ’--[1]-( 5 4 5 7 2| 23)| ’--[5]-( 0 0 3 0 0| 3)-T+--[4]-( 21 23 15 19 27| 105)| ’--[5]-( 6 4 6 3 6| 25)| ’--[2]-( 1 0 0 0 4| 5)-T’--[5]-( 18 22 14 25 24| 103)’--[2]-( 0 5 3 5 1| 14)-TThe resulting context function can be derived from the given R output,it is given below.424.1. SIMULATION STUDIESc(x0−∞) =c2[c2c3c4c5], x0 = c2, x−1 ∈ {c2c3c4c5}, x−2−∞ arbitraryc2c1[c1c2c3c3], x0 = c2, x−1 = c1, x−2 ∈ {c1c2c3c3}, x−3−∞ arbitraryc1c2c5, x0 = c2, x−1 = c2, x−2 = c5, x−3−∞ arbitraryc4[c1c2c3c4], x0 = c4, x−1 ∈ {c1c2c3c4}, x−2−∞ arbitraryc4c5[c1c3c4c5], x0 = c4, x−1 = c5, x−2 ∈ {c1c3c4c5}, x−3−∞ arbitraryc4c5c2, x0 = c4, x−1 = c5, x−2 = c2, x−3−∞ arbitraryc5[c1c3c4c5], x0 = c5, x−1 ∈ {c1c3c4c5}, x−2−∞ arbitraryc5c2, x0 = c5, x−1 = c2, x−2−∞ arbitraryThe context function from the simulated cluster assignments aligns withthe summary output, as 8 contexts are derived, and the maximum historyconsidered is 3 states, meaning an order 3 VLMC.Drawing a VLMC fit in R using draw() from the VLMC package also pro-vides the estimated transition probability vectors for each generated context[20]. In other words, for every Pˆ (xt | cˆ(xt−1, xt−2, . . . )), the values can beobtained directly from the R output. Each line in the output corresponds toan estimated transition probability vector associated with a specific context.For example, the first line, which reads [2]-( 21 15 23 23 14| 96), cor-responds to the context c2[c2c3c4c5], which provides the most likely clusterassignments for the next state, given the current context c2[c2c3c4c5]. Theoutput indicates that given the context c2[c2c3c4c5], the most likely clusterlabel for the upcoming observation is c3 or c4, the transition probabilitiesbelow:Pˆ (Xt = c1 | cˆ(Xt−1 = c2, Xt−2 ∈ {c2c3c4c5}, . . . )) = 2196,Pˆ (Xt = c2 | cˆ(Xt−1 = c2, Xt−2 ∈ {c2c3c4c5}, . . . )) = 1596,Pˆ (Xt = c3 | cˆ(Xt−1 = c2, Xt−2 ∈ {c2c3c4c5}, . . . )) = 2396,Pˆ (Xt = c4 | cˆ(Xt−1 = c2, Xt−2 ∈ {c2c3c4c4}, . . . )) = 2396,Pˆ (Xt = c5 | cˆ(Xt−1 = c2, Xt−2 ∈ {c2c3c4c5}, . . . )) = 1496.This output of the VLMC fit describes the estimated probabilities for thenext state Xt, conditioned on the context history, which includes the rele-vant past observations up to order k. These results can be used to predictfuture states from X , by referencing the estimated transition probabilities434.2. REAL-WORLD APPLICATIONSfor each observed context. This will be extended to our model-based clus-tering results to identify the most likely class for an upcoming calendar yearof weekly wildfire counts for the given context history.4.2 Real-World ApplicationsIn this section, the counts and causes of wildfires are studied with theCNFDB dataset in terms of forecasting and clustering results. The firstapplication will be in Alberta for the south, northeast and northwest regionsin terms of forecast analysis of wildfire causes using VLMC models. The goalis to study the use of VLMC to forecast the causes in unseen data, wherebalanced and unbalanced cases will be explored to validate the extensionto the model-based clustering approach. Additionally, comparing VLMCsto a baseline first order Markov chain in terms of out-of-sample forecastingmetrics is completed. Weekly wildfire counts will also be analyzed in termsof intermediate-term out-of-sample forecasting using both the BATS and theTBATS models for comparison. Finally, model-based clustering of weeklywildfire time series using TBATS models is explored for BC, providing abasis for future work.4.2.1 Predicting Wildfire Causes using VLMCsTo study the cutoff parameter K for VLMC fits, K will be varied whenfitting VLMC models and to evaluate performance by comparing both in-sample and out-of-sample accuracy measures to identify the causes of wild-fires in Alberta.Figure 5.14 provides a visual summary of the distribution of wildfirecauses in Alberta from the CNFDB dataset. Alberta is then divided intothree regions: northeast, northwest, and south for the prediction of wildfirecauses based on Figure 5.14. Each region is divided into a training set (80%)and a test set (20%) to compare a variety of precision metrics for differentwildfire cause distributions. The cutoff parameter K will range from 0.5to 3.5 and metrics such as accuracy, precision, recall, and the F1 score arereported for each class and some averaged metrics. The considered causesof wildfires are human (H) and natural (N), since these account for almostall documented wildfires in Alberta [22]. Furthermore, including unknown(U) causes reduces the out-of-sample metrics, as the class imbalance is quitesevere.For southern Alberta, the number of human-caused wildfires tends tobe more than natural, whereas in northern Alberta natural wildfire causes444.2. REAL-WORLD APPLICATIONStend to occur at a higher rate. Specifically, southern Alberta had 15,957human-caused and 6,162 natural-caused wildfires, northwest Alberta had10,562 human-caused and 10,298 natural-caused wildfires, and northeastAlberta had 4,030 human-caused and 7,586 natural-caused wildfires. Thediscrepancy between the causes of the wildfires in all three regions allows thepredictive capabilities of VLMC to be assessed over a realistic distributionthat can be extended across Canada.For all three regions, increasing the cutoff parameterK led to a reductionin overfitting as the context tree for each VLMC fit began to decrease in size.Specifically, once K increased from 0.5 precision, recall, and F1 score, beganto increase for the test set for both causes, followed by a reduction of thesemetrics for the training set. For all three regions, the optimal value for Kwas around 2-3, which aligns with the recommendation for this parameterin [24].Southern and northeast Alberta exhibited class imbalances in the wildfirecause distribution. Specifically, for southern Alberta, the out-of-sample met-rics for human-caused wildfires were approximately: precision = 0.85, recall= 0.90, and F1 score = 0.87, while the out-of-sample metrics for natural-caused wildfires were approximately: precision = 0.46, recall = 0.35, and F1score = 0.40. The southern region documented more human-caused wildfirescompared to natural-caused ones, as can be seen with the reduced out-of-sample metrics for natural-caused wildfires in this region. In comparison,the out-of-sample metrics for human-caused wildfires in northeast Albertawere approximately: precision = 0.80, recall = 0.75, and F1 score = 0.77,while the out-of-sample metrics for natural-caused wildfires were approxi-mately: precision = 0.75, recall = 0.80, and F1 score = 0.78. In the north-east region of Alberta, the documented number of natural-caused wildfireswas higher than the human-caused, but had slightly less class imbalance.This can be seen with out-of-sample metrics for natural-caused wildfires,as the results improved compared to the southern region. Out-of-samplemetrics for human-caused wildfires in southern Alberta provided good re-sults, while natural-caused metrics were relatively poor. Furthermore, theout-of-sample metrics for human-caused wildfires in northeast Alberta sawa reduction compared to southern Alberta, while the natural-caused metricsshowed a significant improvement compared to the southern region. Consid-ering that both regions of the wildfire cause distributions were unbalanced,predicting from a VLMC fit on unbalanced data can yield different resultsdepending on the specific structure.Taking into account the class imbalance of the wildfire causes for both thesouthern and northeast regions, the out-of-sample metrics can be promising454.2. REAL-WORLD APPLICATIONSand indicate that the model can generalize to complex unseen real noisydata. Specifically, both regions provide varying F1 scores, which also pro-duce macro-averaged F1 scores that demonstrate that VLMCs can workwith class imbalances. Furthermore, the southern region reports recall val-ues that indicate that VLMC can correctly identify human-caused wildfiresconsistently, while low precision values indicate that predicted unknown-caused wildfires are typically human-caused. In contrast, Northeast Albertaproduces slightly lower out-of-sample precision and recall metrics for human-caused wildfires, while major improvements are made for out-of-sample pre-cision and recall metrics for natural-caused wildfires. Taking into accountboth imbalanced regions, northeast Alberta outperforms southern Albertain this study, which may be due to the fact that southern Alberta exhibits alarger class imbalance of human-caused wildfires. Tables 5.7 and Table 5.9provide the complete summary of the southern and northeast regions ofAlberta.The northwest of Alberta exhibited an approximately equal distributionof wildfire causes between human and natural origins. Out-of-sample met-rics were similar to those of the northeast region, while outperforming thesouthern region. Specifically, the out-of-sample metrics for human-causedwildfires were approximately: precision = 0.77, recall = 0.79, and F1 score= 0.78, while the out-of-sample metrics for natural-caused wildfires wereapproximately: precision = 0.73, recall = 0.71, and F1 score = 0.72. Simi-larly to the northeast region of Alberta, the implementation of VLMCs as aforecasting method for wildfire occurrences produces reliable and consistentresults on unseen complex data. The results for this region solidify the ef-fectiveness of VLMCs for data that are not severely imbalanced in terms ofout-of-sample performance, and serves as a viable method to forecast wild-fire causes for different regions in Canada. Table 5.8 provides the detailsdiscussed for the northwest region.4.2.2 VLMC versus First Order Markov Chain: WildfireCause Class ImbalancePredicting from a VLMC model for wildfire causes that can be imbal-anced or balanced provided reliable and varying out-of-sample results, espe-cially for the imbalanced case. This was especially true when the causes ofthe wildfires focused on human and natural causes, which account for themajority of documented cases in Canada. The addition of unknown wild-fires further imbalances the distribution and can reduce the performance ofa VLMC in terms of out-of-sample performance. To assess the viability of464.2. REAL-WORLD APPLICATIONSVLMCs for real data that can exhibit severe class imbalances, a compar-ison will be made to a first order Markov chain in terms of out-of-sampleperformance.The southern, northeast and northwest regions of Alberta will be consid-ered again, with the addition of including unknown causes of wildfires. Dueto the addition of unknown wildfire causes, the out-of-sample metric willdecrease for VLMC forecasting results. The distribution of unknown causesof wildfires is 773 in southern Alberta, 555 in northwest Alberta, and 221in northeast Alberta. In each region of Alberta, unknown wildfire causesaccount for less than 4% of the reported wildfire causes, providing a classthat is represented much less than human and natural causes.The results can be compared to a baseline of a first order Markov chainthat was inspired from modeling of the Fire Weather Index (FWI) separatedinto different classes from [21]. Specifically, the out-of-sample accuracy isreported for each region, which is slightly misleading due to the severe imbal-ance of the causes of the wildfires. The macro-averaged F1 score is reportedas the metric to contrast the VLMCs forecasting results to the first-orderMarkov chain in favor of the accuracy metric. The macro averaged F1 scoreis the average of the F1 scores for each wildfire cause considered for differ-ent VLMC and first order Markov Chain forecasting results for the differentregions. To allow each VLMC fit to grow a large context tree and handlethe class imbalance, the cutoff parameter will be fixed at K = 0.5 for eachregion. Furthermore, a naive random model from the observed trainingproportions of all three causes will be used to compare the out-of-samplemetrics to the VLMC and the first order Markov chain. The naive randommodel results were computed from 10,000 random samples from the sametraining sets as the VLMC and the first order Markov chain.Southern and Northwestern Alberta produced greater macro-averagedF1 scores for the VLMC models in comparison to the first order Markovchains. Specifically, the southern region in Alberta had a noticeable increasein deploying the VLMC forecasting in comparison to the first order MarkovChain for the macro-averaged F1 score. The VLMC in northwest regiononly provided marginal improvements for the macro-averaged F1 score overthe first order Markov chain results. Furthermore, Northeastern Albertaactually saw the first order Markov chain having improved accuracy andmacro-averaged F1 scores in comparison to the VLMC results.The accuracy measures for all the regions in Alberta originating fromthe VLMC fits actually say lower accuracy scores in comparison to the firstorder Markov chain results. This is not alarming, as the wildfire causesdistribution for all three regions were imbalanced. The difference in VLMCs474.2. REAL-WORLD APPLICATIONSand the first order Markov chains in this application revolves around theadvantage of the VLMCs to predict unknown wildfire causes in comparisonto first order Markov chains, where the Markov chains metrics tend to bebiased to human and natural causes. These out-of-sample results illustratethe limitations of both VLMCs and first order Markov Chains when dealingwith severe class imbalances, while also demonstrating that VLMCs canoutperform first order Markov Chains when dealing with class imbalancesin terms of identifying minority classes.The VLMC and the first order Markov chain metrics significantly out-performed the naive random model for all three regions. The naive randommodel assumed independence upon each sample from the trained propor-tions, which indicates the serial dependence approach from the VLMC andfirst order Markov chain a better fit for modeling wildfire causes in Canada.Table 4.2 provides the complete summary of these results.Table 4.2: Comparison of VLMC first order Markov chain to an independentsampling baseline by region for wildfire cause in Alberta. Metrics includeout-of-sample accuracy and macro-averaged F1 score.Region Model Accuracy Macro F1South VLMC 0.7433 0.4758First Order Markov Chain 0.7549 0.4172Independent Sampling 0.5728 0.3169Northwest VLMC 0.6756 0.4671First Order Markov Chain 0.7223 0.4620Independent Sampling 0.4852 0.3467Northeast VLMC 0.7076 0.4802First Order Markov Chain 0.7461 0.5053Independent Sampling 0.4721 0.32094.2.3 Intermediate-Term Forecasting of WildfiresBased on simulation studies completed in British Columbia for the TBATSmodel, forecast accuracy can be compared between the BATS and TBATSmodels. Specifically, British Columbia from 1998-2024 will be studied interms of out-of-sample forecast accuracy between BATS and TBATS modelfits.The training set in British Columbia will include the years 1998-2017,484.2. REAL-WORLD APPLICATIONSwhile the testing set will be 2018-2023. The previous rolling window pro-cedure will be implemented in terms of forecasts for each year in the testset. To evaluate the two models, the out-of-sample mean squared error(MSE) and the mean absolute error (MAE) are computed for each year andaveraged.Figure 4.6: Comparison of BATS (blue) and TBATS (red) out-of-samplemodel performance for weekly wildfire counts in British Columbia (2018-2023). Forecast accuracy is measured using MSE (left panel) and MAE(right panel) from 2018-2023.TBATS outperforms BATS in terms of out-of-sample MSE and MAEperformance on average and for most of the years forecasted. Specifically,the BATS forecasts provide lower MAE values compared to the TBATSforecasts for 2021. However, BATS forecasts have lower MSE for 2018, 2021,and 2023 compared to TBATS forecasts. On average, the out-of-sampleMSE and MAE of TBATS are lower than the predicted weekly wildfirecounts of BATS. This indicates that the additions to the TBATS model fitsprovide more accurate and reliable wildfire forecasts compared to the BATSmodel. Figure 4.6 provides a visual summary, while Table 5.4 records theexact numerical results for more detail.Furthermore, Table 5.3 provides the parameters of the BATS and TBATSmodel for each fit in BC. Repetition of this forecast procedure was also com-pleted to compare the accuracy of the forecast out of the sample for theBATS and TBATS models for some other provinces in Canada. Specifically,similar results to those of British Columbia were found in both Alberta andSaskatchewan, where TBATS outperforms BATS in out-of-sample forecastaccuracy on average. Figures 4.7 and 4.8 provide the results visually, whileTables 5.5 and 5.6 provide the exact numerical summaries for Alberta andSaskatchewan.494.2. REAL-WORLD APPLICATIONSFigure 4.7: Comparison of BATS (blue) and TBATS (red) out-of-samplemodel performance for weekly wildfire counts in Alberta (2018-2023). Fore-cast accuracy is measured using MSE (left panel) and MAE (right panel)from 2018-2023.Figure 4.8: Comparison of BATS (blue) and TBATS (red) out-of-samplemodel performance for weekly wildfire counts in Saskatchewan (2018-2023).Forecast accuracy is measured using MSE (left panel) and MAE (right panel)from 2018-2023.4.2.4 Clustering Wildfire Time SeriesClustering weekly wildfire counts over a calendar year involves extendingthe TBATS model fits. Specifically, Algorithm (1) is executed in which theresulting model fits are passed to Algorithm (2) and the features are ex-tracted and scaled to perform clustering. For a preliminary result, 74 years(1950-2023) of weekly wildfire counts in BC were clustered. Upon fittingeach calendar year to a TBATS model, similar weekly wildfire counts wereassigned to the clusters using a model-based clustering technique, with theEM algorithm for estimation. This process was automated with Mclust()from the mclust package [28]. Once completed, generate a representativetime series for each cluster by averaging the weekly fire counts for the clus-504.2. REAL-WORLD APPLICATIONSter; Fit a TBATS model for each cluster average to summarize typical fireoccurrence patterns for the cluster. The residual diagnostics for the cluster-ing results had some years with residual autocorrelation. These years willbe kept in the analysis to allow the clustering procedure and the VLMCextension to have sufficient data.The 74 years of BC wildfire data were distributed across the clusters asfollows: Cluster 1 had 39 years, Cluster 2 had 9 years, and Cluster 3 had 26years. Once the representative weekly wildfire counts for each cluster wasrecorded, this was then refit as a TBATS model for a model-based clusteringapproach that could be simulated from. The simulation results form threefairly distinct cluster patterns for the weekly wildfire counts, which alignclosely with their average representations. Specifically, cluster 1 providesinformation that the typical pattern has a longer duration of increase inwildfires in BC around early summer, while it forms a steep increase aroundlate summer after a rapid decrease. Cluster 2 provides a somewhat similarpattern to cluster 1, where in this case there is a smaller increase in earlysummer with additional increases around the expected increase in late sum-mer. Finally, cluster 3 provides a pattern with a gradual increase duringthe summer with less significant increases or decreases for the number ofwildfire occurrences in BC, illustrating consistent growth until late summer.All three cluster results illustrate the simulated weekly wildfire counts of thefinal TBATS fit that align with their respective typical pattern; Figures 4.9to 4.11 provides the illustrations for all the cluster results.514.2. REAL-WORLD APPLICATIONSFigure 4.9: TBATS simulated weekly wildfire counts for Cluster 1 in BC,using the average of all years in the cluster (1950–2023).Figure 4.10: TBATS simulated weekly wildfire counts for Cluster 2 in BC,using the average of all years in the cluster (1950–2023).524.2. REAL-WORLD APPLICATIONSFigure 4.11: TBATS simulated weekly wildfire counts for Cluster 3 in BC,using the average of all years in the cluster (1950–2023).Using the clustering results and the associated labels for each year inBC, a VLMC model fit can be deployed to identify the upcoming calendaryear into one of the three clusters. Denote the state space X = {c1, c2, c3}for the cluster labels that align with the methodology and notation in [24].Following the results in BC and using the VLMC package from [20] withK = 2, the following VLMC output was obtained from the clustering labels.[x]-( 38 9 26| 73)+--[1]-( 20 6 12| 38)| ’--[3]-( 5 1 6| 12)| ’--[2]-( 0 0 3| 3)-T’--[3]-( 13 3 10| 26)+--[1]-( 6 1 5| 12)| ’--[2]-( 1 1 0| 2)-T’--[2]-( 3 1 0| 4)>-TUsing the previous output, the context function can for the VLMC fitcan be obtained. The VLMC has an order of 3 and following [24] the corre-sponding context function is:534.2. REAL-WORLD APPLICATIONSc(x0−∞) =c1[c1c2], x0 = c1, x−1 ∈ {c1c2}, x−2−∞ arbitraryc1c3[c1c3], x0 = c1, x−1 = c3, x−2 ∈ {c1c3}, x−3−∞ arbitraryc1c3c2, x0 = c1, x−1 = c3, x−2 = c2, x−3−∞ arbitraryc3[c2c3], x0 = c3, x−1 ∈ {c2c3}, x−2−∞ arbitraryc3c1[c1c3], x0 = c3, x−1 = c1, x−2 ∈ {c1c3}, x−3−∞ arbitraryc3c1c2, x0 = c3, x−1 = c1, x−2 = c2, x−3−∞ arbitraryc3c2, x0 = c3, x−1 = c2, x−2−∞ arbitraryc2, x0 = c2, x−1−∞ arbitraryFrom the fitted VLMC the minimal state space is given with the contextsto identify the most likely cluster for the upcoming year. The distributionof cluster assignments can be summarized in Figure 4.3.Table 4.3: Cluster assignments for BC wildfire occurrence time series from1950-2023.Cluster Yearsc1 1950, 1951, 1954, 1956, 1958, 1960, 1961, 1967, 1972, 1976,1980, 1984, 1987, 1988, 1989, 1991, 1992, 1993, 1994, 1995,1996, 1998, 1999, 2000, 2002, 2003, 2004, 2005, 2007, 2009,2010, 2011, 2012, 2014, 2015, 2016, 2017, 2018, 2019, 2023c2 1952, 1965, 1977, 1979, 1982, 2001, 2006, 2013, 2015c3 1953, 1955, 1957, 1959, 1962, 1963, 1964, 1966, 1968, 1969,1970, 1971, 1973, 1974, 1975, 1978, 1981, 1983, 1985, 1986,1990, 1997, 2008, 2020, 2021, 2022Considering the assignment of clusters, we want to be able to identify thecluster for the upcoming yearly weekly wildfire occurrence. Given the clusterlabels for each year, the most recent context of highest order is c1c3[c1c3],which corresponds to the years 2023, 2023, and 2021. Given this informationand the context tree drawn before from the VLMC package, the most likelycluster label for the upcoming year of weekly wildfire counts is c3. This isfrom the second line in the drawn context tree [3]-( 5 1 6| 12), whichprovides the transition probabilities as:544.2. REAL-WORLD APPLICATIONSPˆ (Xt = c1 | cˆ(Xt−1 = c2, Xt−2 = c3, Xt−3 ∈ {c1c3}, . . . )) = 512,Pˆ (Xt = c2 | cˆ(Xt−1 = c2, Xt−2 = c3, Xt−3 ∈ {c1c3}, . . . )) = 112,Pˆ (Xt = c3 | cˆ(Xt−1 = c2, Xt−2 = c3, Xt−3 ∈ {c1c3}, . . . )) = 12.The model-based clustering approach in BC from 1950-2023 for weeklywildfire occurrences provided preliminary results to help cluster wildfireoccurrence time series. Furthermore, the model-based approach provideda representative TBATS model for each cluster which was extended to aMarkov chain modeling approach to help identify an upcoming wildfire sea-son.55Chapter 5Conclusion5.1 SummaryThis thesis studied the occurrence of wildfires and causes through appli-cations using methods from [8, 24]. Wildfire occurrence forecasting, sim-ulation, and clustering applications were conducted for Alberta, BritishColumbia, and Saskatchewan, using the BATS and TBATS model frame-work in [8]. The motivation for wildfire occurrence applications originatedfrom discussions with members of the British Columbian provincial govern-ment along with the increase of wildfires in Canada [22]. Furthermore, thewildfire cause applications in this thesis originated from a graduate courseproject, where wildfire causes such as human and natural tend to cause themajority destructive wildfires while accounting for the majority [22].The exponential smoothing methods of BATS and TBATS provided astable and promising basis for forecasting, simulating, and clustering theoccurrence of wildfires across various regions of Canada. The frameworkprovided from [8] allowed for likelihood computation, model selection, andmore from time series data for wildfire occurrence. The results of this thesisdemonstrated the effectiveness and stability of TBATS models to simulatewildfire occurrences to training or testing data in British Columbia; In Sec-tion 4, a rolling window out-of-sample simulation study was conducted inBritish Columbia and demonstrated the consistency and precision in captur-ing the overall pattern for weekly wildfire counts. Furthermore, increasingthe number of simulations produced simulated weekly wildfire counts thatconverged to the typical pattern and smoothed out small fluctuations inthe process, while producing confidence intervals that captured the array ofsimulation patterns obtained. This rolling window approach was extendedto compare out-of-sample metrics for BATS and TBATS models on weeklywildfire occurrence data in British Columbia. The results indicated that forAlberta, British Columbia, and Saskatchewan TBATS outperformed BATSon average in reducing MSE and MAE, illustrating that both methods areviable for forecasting wildfire occurrences in Canada.A model-based clustering framework was constructed based on TBATS565.2. FUTURE WORKmodel fits to simulate or forecast wildfire occurrences. The preliminary re-sults demonstrated the typical weekly patterns of wildfire occurrence fora calendar year in British Columbia for each cluster, while presenting asimulated comparison that aligned with the representation of each group.Furthermore, the methodology of VLMCs from [24] was implemented onthe clustering labels for the 74 years of data in BC. The results demon-strate that using our model-based clustering approach in addition to VLMCmodelling, an upcoming wildfire season can be identified from one of the ex-isting clusters. Our results indicated that using data in BC from 1950-2023,the upcoming year of 2024 has the highest likelihood of aligning with theTBATS model of cluster 3.The work on wildfire causes using VLMCs from [24] demonstrated theeffectiveness of the variable-length memory approach compared to finite or-der Markov chains. Additionally, the out-of-sample metrics for both theVLMC and the first order Markov chain in all regions in Alberta highlightthe sequential dependence modelling approach over independent sampling.Furthermore, the VLMC study on wildfire causes provided insight on theperformance of VLMCs on balanced and unbalanced data, indicating themodels effectiveness in handling varying class distributions. Finally, thewildfire cause work leveraged the use of VLMCs to model cluster labels forthe model-based clustering approach.5.2 Future WorkThis thesis formed a basis for a model-based clustering approach to helpidentify an upcoming calendar year as a known class modelled by a TBATSmodel. In the future, the following steps and additions could be done toextend and improve our method to provide more reliable and consistentconclusions.− Study additional regions throughout Canada for the modeling, simu-lation, and forecasting of BATS and TBATS.− Reduce the area of study for any given region. This in turn wouldproduce results that are specific and applicable to a specific region,making an immediate impact.− Add covariates to the analysis for weekly wildfire occurrences. Specif-ically, the use of FWI could enhance the ability of the clustering ap-proach by adding additional information.57Bibliography[1] S. Aghabozorgi, A. Seyed Shirkhorshidi, and T. Ying Wah,Time-series clustering – a decade review, Information Systems, 53(2015), pp. 16–38.[2] T. W. Anderson, The Statistical Analysis of Time Series, John Wiley& Sons, Inc., 2011.[3] M. Auger-Me´the´, K. Newman, D. Cole, F. Empacher,R. Gryba, A. A. King, V. Leos-Barajas, J. Mills Flemming,A. Nielsen, G. Petris, and L. Thomas, A Guide to State–SpaceModeling of Ecological Time Series, Ecological Monographs, 91 (2021),p. e01470.[4] N. Banerjee, A. Morton, and K. Akartunalı, Passenger de-mand forecasting in scheduled transportation, European Journal of Op-erational Research, 286 (2020), pp. 797–810.[5] P. Bu¨hlmann and A. J. Wyner, Variable length Markov chains, TheAnnals of Statistics, 27 (1999), pp. 480 – 513.[6] Canadian Forest Service, Canadian national fire database – agencyfire data. Natural Resources Canada, 2024. Accessed: 2024-07-30.[7] J. D. Cryer and K.-S. Chan, Time series analysis: with applicationsin R, Springer, 2 ed., 2008.[8] A. M. De Livera, R. J. Hyndman, and R. D. Snyder, Forecastingtime series with complex seasonal patterns using exponential smoothing,Journal of the American Statistical Association, 106 (2011), pp. 1513–1527.[9] R. Durrett, Essentials of Stochastic Processes, Springer, 3 ed., 2012.[10] H. Economics, The full community costs of wildfire, Headwaters Eco-nomics: Bozeman, MT, USA, (2018), p. 44.58[11] E. S. Gardner and E. Mckenzie, Forecasting trends in time series,Management Science, 31 (1985), pp. 1237–1246.[12] S. Ghassempour, F. Girosi, and A. Maeder, Clustering multivari-ate time series using hidden markov models, International Journal ofEnvironmental Research and Public Health, 11 (2014), pp. 2741–2763.[13] L. Han, Time Series Course Notes. Unpublished manuscript, January2024.[14] C. C. Hanes, X. Wang, P. Jain, M.-A. Parisien, J. M. Little,and M. D. Flannigan, Fire-regime changes in Canada over the lasthalf century, Canadian Journal of Forest Research, 49 (2019), pp. 256–269.[15] R. Hyndman, G. Athanasopoulos, C. Bergmeir, G. Caceres,L. Chhay, M. O’Hara-Wild, F. Petropoulos, S. Razbash,E. Wang, and F. Yasmeen, forecast: Forecasting functions for timeseries and linear models, 2024. R package version 8.23.0.[16] R. J. Hyndman and G. Athanasopoulos, Forecasting: Principlesand Practice, OTexts, Melbourne, Australia, 2nd ed., 2018. https://otexts.com/fpp2/.[17] R. J. Hyndman, A. B. Koehler, R. D. Snyder, and S. Grose,A state space framework for automatic forecasting using exponentialsmoothing methods, International Journal of Forecasting, 18 (2002),pp. 439–454.[18] M. Khodarahmi and V. Maihami, A review on Kalman filter models,Archives of Computational Methods in Engineering, 30 (2023), pp. 727–747.[19] G. Kitagawa, Introduction to Time Series Modeling, CRC Press, 2010.[20] M. Maechler, VLMC: Variable Length Markov Chains (’VLMC’)Models, 2024. R package version 1.4-4.[21] D. L. Martell, A Markov chain model of day to day changes in thecanadian forest fire weather index, International Journal of WildlandFire, 9 (1999), pp. 265–273.[22] T. McGee, B. McFarlane, and C. Tymstra, Chapter 3 - wildfire:A canadian perspective, in Wildfire Hazards, Risks and Disasters, J. F.Shroder and D. Paton, eds., Elsevier, Oxford, 2015, pp. 35–58.59[23] Ministers Responsible for Emergency Management,An emergency management framework for Canada, 2017.https://www.publicsafety.gc.ca/cnt/rsrcs/pblctns/2017-mrgnc-mngmnt-frmwrk/2017-mrgnc-mngmnt-frmwrk-en.pdf.[24] M. Ma¨chler and P. Bu¨hlmann, Variable length Markov chains:Methodology, computing, and software, Journal of Computational andGraphical Statistics, 13 (2004), pp. 435–455.[25] F. Orderud, Comparison of Kalman filter estimation approaches forstate space models with nonlinear measurements, in Proc. of Scandina-vian Conference on Simulation and Modeling, 2005, pp. 1–8.[26] S. S. Rangapuram, M. W. Seeger, J. Gasthaus, L. Stella,Y. Wang, and T. Januschowski, Deep state space models for timeseries forecasting, in Advances in Neural Information Processing Sys-tems, S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi, and R. Garnett, eds., vol. 31, Curran Associates, Inc., 2018.[27] L. Scrucca, M. Fop, T. B. Murphy, and A. E. Raftery, mclust5: clustering, classification and density estimation using gaussian finitemixture models, The R journal, 8 (2016), p. 289.[28] L. Scrucca, C. Fraley, T. B. Murphy, and A. E. Raftery,Model-Based Clustering, Classification, and Density Estimation Usingmclust in R, Chapman and Hall/CRC, 2023.[29] R. H. Shumway and D. S. Stoffer, Time Series Analysis and ItsApplications: with R Examples, Springer, 4 ed., 2017.[30] J. W. Taylor, Short-term electricity demand forecasting using doubleseasonal exponential smoothing, Journal of the Operational ResearchSociety, 54 (2003), pp. 799–805.[31] H. Teichgraeber and A. R. Brandt, Clustering methods to findrepresentative periods for the optimization of energy systems: An initialframework and comparison, Applied Energy, 239 (2019), pp. 1283–1293.[32] R. S. Tsay, Time series and forecasting: Brief history and futureresearch, Journal of the American Statistical Association, 95 (2000),pp. 638–643.60[33] C. Tymstra, B. J. Stocks, X. Cai, and M. D. Flannigan, Wild-fire management in Canada: Review, challenges and opportunities,Progress in Disaster Science, 5 (2020), p. 100045.[34] T. Ullrich, On the autoregressive time series model using real andcomplex analysis, Forecasting, 3 (2021), pp. 716–728.[35] R. Umatani, T. Imai, K. Kawamoto, and S. Kunimasa, Timeseries clustering with an em algorithm for mixtures of linear gaussianstate space models, Pattern Recognition, 138 (2023), p. 109375.[36] B. Venkataramana Kini and C. Chandra Sekhar, Bayesian mix-ture of ar models for time series clustering, Pattern Analysis and Ap-plications, 16 (2013), pp. 179–200.[37] N. Vlassis and A. Likas, A greedy em algorithm for gaussian mixturelearning, Neural processing letters, 15 (2002), pp. 77–87.[38] G. Welch, G. Bishop, et al., An introduction to the Kalman filter,(1995).[39] M. West and J. Harrison, Bayesian forecasting and dynamic mod-els, Springer, 2 ed., 1997.[40] P. R. Winters, Forecasting sales by exponentially weighted movingaverages, Management Science, 6 (1960), pp. 324–342.[41] B. M. Wotton, C. A. Nock, and M. D. Flannigan, Forest fire oc-currence and climate change in Canada, International Journal of Wild-land Fire, 19 (2010), pp. 253–271.[42] Y. Xiong and D.-Y. Yeung, Time series clustering with ARMAmixtures, Pattern Recognition, 37 (2004), pp. 1675–1689.61AppendicesAppendix A: TablesTable 5.1: Notable variables and descriptions from the Canadian NationalFire Database (CNFDB) dataset.Variable DescriptionFIRE ID Fire agency identificationSRC AGENCY Agency fire data was obtained (province, territory, or parks)NAT PARK Name of national park affected (if applicable)LATITUDE Latitude of the fire start locationLONGITUDE Longitude of the fire start locationYEAR Year for the fire provided from agenciesMONTH Month for the fire provided from agenciesDAY Day for the fire provided from agenciesREP DATE Date for the fire from agency (YEAR-MONTH-DAY)OUT DATE Date the fire was reported as outSIZE HA Size of fire in hectares reportedCAUSE Reported cause of fireRESPONSE Type of response initiated by the agencyFIRE TYPE Agency reported fire type62Table 5.2: Components of BATS and TBATS model fits. TBATS includesadditional components for the seasonal components.Component Descriptionmodel Fitted modely Input seriesfitted.values Fitted valueserrors Model residualslambda Box–Cox parameterdamping.parameter Damping parameteralpha Level smoothing parameterbeta Trend smoothing parametergamma.one.values First seasonal smoothing parameter (TBATS)gamma.two.values Second seasonal smoothing parameter (TBATS)gamma.values Seasonal smoothing parameters (BATS)ar.coefficients Autoregressive coefficientsma.coefficients Moving average coefficientsseasonal.periods Seasonal periodsk.vector Fourier harmonics per seasonal periodvariance Variance estimatelikelihood Log-likelihood of the modelAIC Akaike Information Criterion63Table 5.3: Model summary of fitted BATS and TBATS models for weeklywildfire counts in British Columbia (2018-2023).Year BATS Model Fit TBATS Model Fit2018 BATS{1, {1, 4}, 0.946, 52} TBATS{1, {2, 2},−, {52.18, 6}}2019 BATS{1, {1, 4}, 0.977, 52} TBATS{1, {2, 2},−, {52.18, 6}}2020 BATS{1, {1, 4}, 0.987, 52} TBATS{1, {0, 3},−, {52.18, 6}}2021 BATS{1, {1, 4}, 0.937, 52} TBATS{1, {0, 3},−, {52.18, 6}}2022 BATS{1, {1, 4}, 0.966, 52} TBATS{1, {0, 3},−, {52.18, 7}}2023 BATS{1, {1, 4}, 0.930, 52} TBATS{1, {0, 3},−, {52.18, 6}}Table 5.4: Comparison of BATS and TBATS out-of-sample model perfor-mance for weekly wildfire counts in British Columbia (2018-2023). Forecastaccuracy is measured using MSE and MAE for each year in the test set.Year BATS MSE BATS MAE TBATS MSE TBATS MAE2018 5376.54 27.78 5854.38 24.732019 1016.94 15.89 554.27 13.072020 1344.31 18.76 488.40 12.102021 2881.87 19.20 3367.10 21.012022 3909.06 29.63 3139.55 22.582023 2582.90 23.00 3282.93 22.21Average 2851.94 22.38 2781.11 19.2864Table 5.5: Comparison of BATS and TBATS out-of-sample model perfor-mance for weekly wildfire counts in Alberta (2018-2023). Forecast accuracyis measured using MSE and MAE for each year in the test set.Year BATS MSE BATS MAE TBATS MSE TBATS MAE2018 822.24 14.60 1049.8 15.742019 412.84 11.87 327.08 12.032020 239.70 10.43 193.75 9.232021 917.63 15.72 771.90 14.642022 1476.43 23.76 520.26 13.812023 466.85 13.91 560.98 12.97Average 722.62 15.05 570.64 13.07Table 5.6: Comparison of BATS and TBATS out-of-sample model perfor-mance for weekly wildfire counts in Saskatchewan (2018-2023). Forecastaccuracy is measured using MSE and MAE for each year in the test set.Year BATS MSE BATS MAE TBATS MSE TBATS MAE2018 112.87 5.65 142.43 6.152019 36.47 3.16 35.72 3.282020 40.59 4.17 25.73 3.582021 488.89 8.55 444.18 7.532022 120.87 6.77 86.15 5.122023 80.62 4.95 95.84 4.98Average 146.72 5.54 138.34 5.1165Table 5.7: Out-of-sample VLMC model metrics for southern Alberta by cause. Wildfire cause distribution for theSouth region: H: 15,957; N: 6,162.K ClassTrain TestAccuracy Precision Recall F1 Macro F1 Accuracy Precision Recall F1 Macro F10.5 H 0.8970 0.8958 0.9654 0.9293 0.8697 0.7608 0.8347 0.8747 0.8542 0.59400.5 N 0.8970 0.9006 0.7361 0.8101 0.8697 0.7608 0.3738 0.3015 0.3338 0.59401.0 H 0.8397 0.8582 0.9241 0.8899 0.7974 0.7764 0.8399 0.8908 0.8646 0.61181.0 N 0.8397 0.7824 0.6412 0.7048 0.7974 0.7764 0.4172 0.3151 0.3590 0.61181.5 H 0.8222 0.8456 0.9134 0.8782 0.7747 0.7863 0.8433 0.9007 0.8711 0.62411.5 N 0.8222 0.7492 0.6079 0.6712 0.7747 0.7863 0.4483 0.3254 0.3771 0.62412.0 H 0.8124 0.8423 0.9014 0.8708 0.7642 0.7897 0.8474 0.8995 0.8727 0.63442.0 N 0.8124 0.7225 0.6034 0.6576 0.7642 0.7897 0.4614 0.3470 0.3961 0.63442.5 H 0.8042 0.8391 0.8918 0.8646 0.7553 0.7938 0.8541 0.8956 0.8744 0.64972.5 N 0.8042 0.7018 0.5983 0.6459 0.7553 0.7938 0.4767 0.3834 0.4250 0.64973.0 H 0.8008 0.8383 0.8872 0.8621 0.7520 0.7940 0.8534 0.8970 0.8747 0.64853.0 N 0.8008 0.6928 0.5979 0.6419 0.7520 0.7940 0.4771 0.3788 0.4223 0.64853.5 H 0.8001 0.8382 0.8861 0.8615 0.7513 0.7954 0.8540 0.8981 0.8755 0.65043.5 N 0.8001 0.6908 0.5981 0.6411 0.7513 0.7954 0.4813 0.3811 0.4254 0.650466Table 5.8: Out-of-sample VLMC model metrics for northwest Alberta by cause. Wildfire cause distribution forthe Northwest region: H: 10,562; N: 10,298.K ClassTrain TestAccuracy Precision Recall F1 Macro F1 Accuracy Precision Recall F1 Macro F10.5 H 0.8990 0.8631 0.9455 0.9024 0.8989 0.7111 0.7242 0.7763 0.7493 0.70420.5 N 0.8990 0.9414 0.8537 0.8954 0.8989 0.7111 0.6918 0.6294 0.6591 0.70421.0 H 0.8410 0.8295 0.8536 0.8414 0.8410 0.7327 0.7491 0.7810 0.7647 0.72761.0 N 0.8410 0.8529 0.8288 0.8407 0.8410 0.7327 0.7100 0.6721 0.6905 0.72761.5 H 0.8271 0.8185 0.8351 0.8267 0.8271 0.7451 0.7642 0.7836 0.7738 0.74101.5 N 0.8271 0.8358 0.8193 0.8275 0.8271 0.7451 0.7199 0.6969 0.7082 0.74102.0 H 0.8190 0.8113 0.8254 0.8183 0.8190 0.7535 0.7730 0.7884 0.7806 0.74972.0 N 0.8190 0.8267 0.8127 0.8196 0.8190 0.7535 0.7280 0.7099 0.7188 0.74972.5 H 0.8131 0.8022 0.8250 0.8134 0.8131 0.7593 0.7767 0.7961 0.7863 0.75542.5 N 0.8131 0.8244 0.8014 0.8127 0.8131 0.7593 0.7362 0.7131 0.7245 0.75543.0 H 0.8100 0.8078 0.8074 0.8076 0.8100 0.7583 0.7830 0.7823 0.7826 0.75533.0 N 0.8100 0.8122 0.8126 0.8124 0.8100 0.7583 0.7275 0.7283 0.7279 0.75533.5 H 0.8086 0.8055 0.8074 0.8064 0.8086 0.7593 0.7829 0.7849 0.7839 0.75613.5 N 0.8086 0.8117 0.8097 0.8107 0.8086 0.7593 0.7295 0.7272 0.7283 0.756167Table 5.9: Out-of-sample VLMC model metrics for northeast Alberta by cause. Wildfire cause distribution forthe Northeast region: H: 4,031; N: 7,585.K ClassTrain TestAccuracy Precision Recall F1 Macro F1 Accuracy Precision Recall F1 Macro F10.5 H 0.9052 0.8181 0.8879 0.8519 0.8774 0.7291 0.7394 0.7238 0.7315 0.72540.5 N 0.9052 0.9486 0.9128 0.9304 0.8774 0.7291 0.7188 0.7346 0.7266 0.72541.0 H 0.8523 0.7903 0.7045 0.7450 0.7910 0.7593 0.7896 0.7196 0.7527 0.74681.0 N 0.8523 0.8755 0.9175 0.8961 0.7910 0.7593 0.7329 0.8005 0.7653 0.74681.5 H 0.8393 0.7679 0.6813 0.7212 0.7742 0.7644 0.7946 0.7255 0.7584 0.74191.5 N 0.8393 0.8660 0.9091 0.8870 0.7742 0.7644 0.7381 0.8049 0.7698 0.74192.0 H 0.8297 0.7434 0.6778 0.7097 0.7636 0.7730 0.7946 0.7483 0.7705 0.76332.0 N 0.8297 0.8631 0.8967 0.8796 0.7636 0.7730 0.7531 0.7988 0.7754 0.76332.5 H 0.8269 0.7374 0.6750 0.7050 0.7615 0.7752 0.7977 0.7492 0.7717 0.76792.5 N 0.8269 0.8617 0.8939 0.8773 0.7615 0.7752 0.7545 0.8023 0.7777 0.76793.0 H 0.8207 0.7721 0.5882 0.6652 0.7488 0.7817 0.8322 0.7162 0.7687 0.75533.0 N 0.8207 0.8355 0.9234 0.8772 0.7488 0.7817 0.7421 0.8497 0.7893 0.75533.5 H 0.8207 0.7721 0.5882 0.6652 0.7488 0.7817 0.8322 0.7162 0.7687 0.75533.5 N 0.8207 0.8355 0.9234 0.8772 0.7488 0.7817 0.7421 0.8497 0.7893 0.755368Appendix B: FiguresFigure 5.1: Simulated vs observed wildfire counts in BC for 2018 (Year 1),using TBATS (5 iterations). Trained up to 2017, 95% CI shaded.69Figure 5.2: Simulated vs observed wildfire counts in BC for 2019 (Year 2),using TBATS (5 iterations). Trained up to 20218, 95% CI shaded.Figure 5.3: Simulated vs observed wildfire counts in BC for 2020 (Year 3),using TBATS (5 iterations). Trained up to 2019, 95% CI shaded.70Figure 5.4: Simulated vs observed wildfire counts in BC for 2021 (Year 4),using TBATS (5 iterations). Trained up to 2020, 95% CI shaded.Figure 5.5: Simulated vs observed wildfire counts in BC for 2022 (Year 5),using TBATS (5 iterations). Trained up to 2021, 95% CI shaded.71Figure 5.6: Simulated vs observed wildfire counts in BC for 2023 (Year 6),using TBATS (5 iterations). Trained up to 2022, 95% CI shaded.Figure 5.7: Simulated vs observed wildfire counts in BC for 2018 (Year 1),using TBATS (50 iterations). Trained up to 2017, 95% CI shaded.72Figure 5.8: Simulated vs observed wildfire counts in BC for 2019 (Year 2),using TBATS (50 iterations). Trained up to 2018, 95% CI shaded.Figure 5.9: Simulated vs observed wildfire counts in BC for 2020 (Year 3),using TBATS (50 iterations). Trained up to 2019, 95% CI shaded.73Figure 5.10: Simulated vs observed wildfire counts in BC for 2021 (Year 4),using TBATS (50 iterations). Trained up to 2020, 95% CI shaded.Figure 5.11: Simulated vs observed wildfire counts in BC for 2022 (Year 5),using TBATS (50 iterations). Trained up to 2021, 95% CI shaded.74Figure 5.12: Simulated vs observed wildfire counts in BC for 2023 (Year 6),using TBATS (50 iterations). Trained up to 2022, 95% CI shaded.Figure 5.13: Weekly wildfire occurrences in BC from 1998-2024. The datawas sourced from the CNFDB dataset.75Figure 5.14: Documented wildfire causes {H,N,U} from 1959–2023 in Al-berta, from the CNFDB dataset, using relative information.76Appendix C: AR(1) Model DerivationsUsing equation (2.8), along with linearity, independence, and recursion,the following derivations obtain the mean and variance of a zero-mean AR(1)process, which align with the conclusion in [13].C.1 Expected Value of an AR(1) ProcessE[zt] = E[ϕzt−1 + εt]= ϕE[zt−1] + E[εt]= ϕ E[ϕzt−2 + εt−1]= ϕ2 E[zt−2] + ϕE[εt−1]=...= ϕn E[zt−n] +n−1∑i=0ϕiE[εt−i](5.1)Recursively expand in n steps until reaching the initial time point t − n,which for convenience is denoted as z0.Assuming that the white noise εt has zero mean, i.e., E[εt] = 0 for all t,this simplifies toE[zt] = ϕnE[z0] = 0 (5.2)C.2 Variance of an AR(1) ProcessVar(zt) = Var(ϕzt−1 + εt)= ϕ2 Var(zt−1) + σ2ε= ϕ2(ϕ2Var(zt−2) + σ2ε)+ σ2ε= ϕ4Var(zt−2) + ϕ2σ2ε + σ2ε=...= ϕ2tVar(z0) + σ2εt−1∑i=0ϕ2i.(5.3)Assuming Var(z0) = σ2z and stationarity, as t→∞ and |ϕ| < 1, the varianceconverges to77Var(zt) = σ2z =σ2ε1− ϕ2 . (5.4)C.3 Autocovariance of an AR(1) Processγk = E[ztzt−k]= E[(ϕzt−1 + εt)zt−k]= ϕ E[zt−1zt−k] + E[εt]E[zt−k]= ϕ γk−1 + 0= ϕ E[(ϕzt−2 + εt−1)zt−k]= ϕ2 γk−2=...= ϕkγ0= ϕkσ2z .(5.5)C.4 Autocorrelation of an AR(1) Processρk =Cov(zt, zt−k)Var(zt)=γkσ2z=ϕkσ2zσ2z= ϕk. (5.6)78Appendix D: BATS and TBATSThe logarithmic likelihood for the observed series is given by (2.29).First, recall the density of the observed series given by (2.28) provided in [8].The authors in [8] concentrate on σ2 with its maximum likelihood estimate.σˆ2 = n−1n∑t=1ε2t . (5.7)To determine the log-likelihood we take the logarithm of (2.28) and sub-stitute (5.7). Taking the logarithm of the density of the observed seriesyieldslog(p(yt | x0, ϑ, σ2)) = log(1(2πσ2)n/2exp(− 12σ2n∑t=1ε2t) n∏t=1yω−1t)= log(1(2πσ2)n/2)+ log(exp(− 12σ2n∑t=1ε2t))+ log( n∏t=1yω−1t)= −n2log(2πσ2)− 12σ2n∑t=1ε2t + log(yω−11 ) + · · ·+ log(yω−1n )= −n2log(2πσ2)− 12σ2n∑t=1ε2t + (ω − 1) log(y1) + · · ·+ (ω − 1) log(yn)= −n2log(2πσ2)− 12σ2n∑t=1ε2t + (ω − 1)n∑t=1log(yt)(5.8)Now substitute (5.7) into (5.8) and multiply the result by −2. Thisyields the following log-likelihood used for BATS and TBATS models.log(p(yt | x0, ϑ)) = L∗(x0, ϑ)= −n2log(2π(n−1n∑t=1ε2t))− 12n−1∑nt=1 ε2tn∑t=1ε2t + (ω − 1)n∑t=1log(yt)= −n2log(2π)− n2log(n−1n∑t=1ε2t)− n2+ (ω − 1)n∑t=1log(yt)= −n2log(2π)− n2log( n∑t=1ε2t)+n2log(n)− n2+ (ω − 1)n∑t=1log(yt)79=:n log(2π) + n log( n∑t=1ε2t)−:n log(n) +>n− 2(ω − 1)n∑t=1log(yt)= n log( n∑t=1ε2t)− 2(ω − 1)n∑t=1log(yt). (5.9)80
Thesis/Dissertation
eng
Mathematics
Vancouver : University of British Columbia Library
Attribution-NonCommercial-NoDerivatives 4.0 International
Graduate
Braun, John, 1963-
Han, Lengyi
Clustering wildfire occurrence time series
Text