The Impact of Climatic Factors on Respiratory Pharmaceutical Demand: A Comparison of Forecasting Models for Greece
Viviana Schisa [∗] and Matteo Farn`e
Department of Statistical Sciences, Alma Mater Studiorum - University of Bologna, Bologna, Italy
Abstract
Climate change is increasingly recognized as a driver of health-related outcomes, yet its impact on pharmaceutical demand remains largely understudied. As environmental conditions evolve and extreme weather events intensify, anticipatingtheir influence on medical needs is essential for designing resilient healthcare systems. This study examinesthe relationship between climate variability and the weekly demand for respiratory prescription pharmaceuticals in Greece, based on a dataset spanning seven and a half years. Granger causality spectra are employed to explore potential causal relationships. Following variable selection, four forecasting models are implemented: Prophet, a Vector Autoregressive model with exogenous variables (VARX), Random Forest with Moving Block Bootstrap (MBB-RF), and Long Short-Term Memory (LSTM) networks. The MBB-RF model achievesthe best performance in relative error metrics while providing robust insightsthrough variable importance rankings. The LSTM model outperforms most metrics, highlighting its ability to capture nonlinear dependencies. The VARX model, which includes Prophet-based exogenous inputs, balances interpretability and accuracy, although it is slightly less competitive in overall predictive performance.
These findings underscore the added value of climate-sensitive variables in modeling pharmaceutical demand and provide a data-driven foundation for adaptive strategies in healthcare planning under changing environmental conditions. 1 Introduction
Climate change represents one of the foremost challenges of recent decades, with profound impacts on environmental systems, public health, and the global economy. This phenomenon is projected to exacerbate global inequality and public health risks, pushing more than 100 million into poverty and limiting accessto vital resources such as food, clean water, and healthcare. These deprivations are expected to raise both mortality and morbidity rates. Accordingto WHO projections, climate change could be responsible for 250,000 additional deaths annually between 2030 and 2050, primarily through direct healththreats. Furthermore, climate change could increase global mortality risk by approximately 3.2% of global GDP by 2100 under a high emissions scenario accordingto projections by Carleton et al.. ∗ Corresponding author: viviana.schisa2@unibo.it 1 Damages are expressed in terms of the Value of a Statistical Life (VSL), an economic measure of the value placed on reducingthe risk of death. The VSL is derived from the willingnessto pay for marginal reductions in mortality risk. It is used in cost-benefit analysesto assessthe economic impacts of policiesthat affect healthand safety. Beyond macroeconomic estimates, climate change is most notably expected to reshape individual healthtrajectories and treatment needs, altering disease incidence and drug demand, particularly in the context of rising chronic conditions and increased exposure to environmental stressors. Due to global warming, the Northern Hemisphere is expected to face not only hotter summers but also more extreme winter conditions. Although average winter temperatures may rise, climate projections indicate an increased likelihood of severe winter storms. These events often bring strong windsthat amplify cold exposure, particularly among vulnerable populations, thereby increasingthe risk of cold-related conditions such as hypothermia and frostbite. In addition, prolonged cold stress can impair immune defenses, making individuals more susceptible to respiratory infections. Risingtemperatures, related increased air pollution and extended pollen seasons intensify respiratory and cardiovascular conditions. In parallel, water and soil quality degradation contribute to arise in diarrheal diseases, especially in vulnerable populations. Furthermore, the spread of vector-borne diseases is extending into non-subtropical regions, and rising anxiety, depression, and stress (related to extreme weather events, displacement, and environmental degradation) are increasingly associated with climate-related stressors and environmental disruptions. These diverse and compounding health stressors are expected to translate into increased pharmaceutical needs across populations. The multifaceted impacts of climate on both physical and mental health are expected to increase the consumption of medications for bothtargeted treatments (e.g., antiprotozoals, psychotropics) and general symptom management (e.g., analgesics). Yet, the relationship between environmental stressors and pharmaceutical demand remains underexplored, even though pharmaceutical consumption could serve as a valuable proxy for climate-related health impacts. Climate change generally refersto long-term transformations in climate systems, but its health impacts can materialize through more immediate and localized events. This study focuses on short- to medium-term fluctuations in environmental conditions, which can be considered atangible expression of climate change. In this sense, increased medication use may reflect heightened exposure to climate-related stressors, providing a quantifiable signal of population-level health impacts. Understandingthese dynamics is crucial to anticipate healthcare needs and support pharmaceutical system resilience. Forecasting models can inform preparedness strategies, optimize resource allocation, and improve supply chain performance, particularly during climate-induced demand surges. Integrated forecasting models incorporating climate, epidemiological, and pharmaceutical data are essential to support policy responses and ensure equitable accessto care. These complex interactions highlight the need for integrated surveillance systems, adaptation strategies, and mitigation effortsthat can deliver both environmental and immediate health benefits. Accurate forecasting enables a proactive strategy to identify vulnerabilities and enhance the robustness of pharmaceutical distribution, e.g., through stockpiling or diversification of high-demand pharmaceuticals. In this context, the present study adopts a quantitative, data-driven approachto forecast respiratory medication demand under climate variability, withthe aim ofinforming both public health policy and pharmaceutical supply chain design. Supply chain resilience can be enhanced by anticipating consumption trends in response to climate change and optimizing resource allocation and asset management, particularly during crises such asthe COVID-19 pandemic. Nonetheless, assessingthe impact of climate change on medical needs remains complex, as healthcare systems respond through a mix of short-term interventions and long-term policy adaptations—mechanismsthat are often multifactorial, delayed, and mediated by socioeconomic conditions, infrastructure constraints, and limited accessto granular health data. Against this backdrop, this study seeksto quantify the impact of climate change on pharmaceutical demand, with a specific focus on respiratory treatments. Amongthe various health domains, respiratory diseases are especially climate-sensitive. Accordingto D’amato et al., global warming accentuatestemperature variability, leadingto higher concentrations of pollutants and extended pollen seasons, which increase allergen exposure and result in more frequent and severe asthma symptoms. Similarly, a survey conducted among physicians from the National Medical Association revealed that patients with COPD and asthma have experienced increased symptoms due to weather changes.
In fact, while direct health impacts, such as mortality, increase due to climate change, adaptation measures, including improved healthcare access and better pharmaceutical supply chains, can mitigate these effects. This study focuses on Greece as a case study, leveraging weekly data over a seven-and-a-halfyear period to explore the relationship between climate variability and respiratory pharmaceutical demand. First, the potential causal links between environmental variables and pharmaceutical consumption are investigated using Granger causality spectra, which allow usto examine frequencyspecific interactions beyond standard time-domain approaches. Based on these insights, multiple forecasting models are then developed and compared —including Prophet, a Vector Autoregressive model with exogenous variables (VARX), Random Forest with Moving Block Bootstrap (MBBRF), and Long Short-Term Memory (LSTM) neural networks—to evaluate the predictive power of climate variables and identify the most effective strategies for anticipating demand fluctuations. The results aim to support public health preparedness and pharmaceutical supply chain resilience under changing climatic conditions. 2 Data
Weekly pharmaceutical sales and climate data for Greece were used to develop time series models and machine learning algorithms, assessing climate-health relationships and improving future healthcare predictions. Alira Health provided data on sales of respiratory pharmaceuticals in Greece from January 4, 2016, to June 26, 2023. Pharmaceutical sales data were obtained through Sell-Out Data Acquisition, enabled by software integration withthe Alira Health panel of 1,200
2 Alira Health, `[https://www.alirahealth.com/](https://www.alirahealth.com/)` pharmacies evenly distributed across Greece. This procedure ensures data quality and reliability while removing all personally identifiable information in full compliance withthe General Data Protection Regulation (GDPR). The statistical expansion was performed weekly, with geographic granularity acrosstwenty distinct territories. It relied on a dynamic panel of pharmacies and incorporated socio-demographic factorsto enhance the accuracy of the estimates. Respiratory drugs classified at the ATC1 level—encompassingthe full spectrum of medicationstargetingthe respiratory system—were used for model estimation and training. This choice was made to account for potential overlaps in therapeutic applications across different respiratory conditions. Following expert advice from Alira, who are highly familiar withthe dataset, only prescription (Rx) medications were considered, asthey are subject to stricter regulatory control and monitoring, thus offering more reliable data for analysis. This aggregated approach does not preclude the possibility of future investigations into specific therapeutic classes or active ingredients. Throughout the manuscript, the weekly time series of respiratory prescription pharmaceuticals sold in Greece will be referred to as drug demand to streamline the discussion. The statistical properties of the series were evaluated usingthe Augmented Dickey–Fuller test with lag order selected by the Bayesian Information Criterion (BIC). The test rejected the null hypothesis of a unit root under the specification with drift, indicatingthat the process is stationary around anon-zero mean. Further analysis of the univariate time series of respiratory drug demand is reported in the Supplementary material, Section S3. The drug demand time series exhibits structural breaks correspondingto the onset and the resolution of the COVID-19 pandemic emergency. These breakpoints were statistically identified usingthe Bai and Perron procedure, which enables consistent estimation of multiple unknown breakpoints within a linear framework. Althoughthe demand for ventilators and face masks increased substantially duringthis period, posing major challenges for the pharmaceutical and medical supply industries, public demand for respiratory drugs decreased significantly.
This decline was largely due to mitigation measures, such as lockdowns and maskwearing, which reduced exposure to environmental, infectious, and allergenic triggers of respiratory diseases, ultimately resulting in structural shifts in pharmaceutical consumption patterns. The climate data, primarily sourced from the ERA5 reanalysis dataset, are available at hourly resolution. Data were collected at the geographical coordinates of the most populous city in each of the twenty pharmaceutical regions of Greece. A regional map showingthe locations used for climate data extraction is provided in the Supplementary material, Section S1. The extracted variables include the 2-meter air temperature ( Temp ), provided in Kelvin and converted to Celsius; the eastward and northward components of the wind at 10 meters above the surface of the Earth; total precipitation ( P ), representing accumulated liquid and frozen water, measured in meters and converted to millimeters; specific humidity, defined asthe mass of water vapour per kilogram of moist air; and proportion of cloud cover.
Wind speed was included because of its influence on pollutant and allergen transport and dispersion. It is derived for each hour hand each region r as
WindSpeed r,h = ~~�~~ u 10 ,r,h [+] [ v] 10 ,r,h [.]
The data were aggregated at the daily and weekly levels, and subsequently at the national level usingthe mean, except for precipitation, which was aggregated usingthe sum to account for the
3 The Anatomical Therapeutic Chemical (ATC) classification system is maintained by the WHO Collaborating Centre: `[https://www.whocc.no/atcdddindex/](https://www.whocc.no/atcdddindex/)` . Level 1 refersto the anatomical main group, e.g., R for the respiratory system. 4 `` Weekly Rx Pharmaceutical Demand
300000
250000
200000
150000
2021−09−20
~~2020−04−06~~
100000
Week start date
Figure 1: Weekly time series of prescription drug demand for respiratory system drugs in Greece, from January 4, 2016, to June 26, 2023. The red dashed lines indicate the mean levels within each segment defined by structural breaks, while the red week-start dates mark the breakpointsthemselves, identified on April 6, 2020, and September 20, 2021. These structural changes were detected usingthe Bai–Perron multiple breakpoint algorithm, withthe optimal number of breaks selected accordingto the BIC. cumulative nature of rainfall over time. To capture short-term climate variability, the number of wet days, extreme rainfall, and temperature standard deviation were derived for each week wand each region r as additional variables from ERA5 data as
= WetDays r,w
= ExtremeRainfall r,w
= TemperatureSD r,w
7 � H
d =1
~~�~~ ~~�~~ �
7
�
7 �( Temp r,d − Temp r,w )
d =1
7 � H × Pr,d
d =1
where Histhe Heaviside step function, P 99 . 9 ,r isthe 99.9th percentile of historical daily precipitation for region r, 1 mm isthe threshold used to define wet days, and Temp r,w isthe regional weekly average temperature. To spatially aggregate these variables across regions, wet days and temperature standard deviation were averaged, asthey reflect relative intensity or frequency, while extreme rainfall was summed to represent the total burden of extreme precipitation events acrossthe country. In addition to the variables mentioned above, the Fire Weather Index (FWI), developed by the Canadian Forestry Service, is sourced from the Copernicus Climate Data Store. The FWI combinesthe Initial Spread Index (ISI), which reflectsthe potential for fire spread, and the Buildup Index (BUI), which indicatesthe amount of fuel available for combustion, thus providing an integrated measure of forest fire conditions, including boththe likelihood of ignition and the potential intensity and spread of fire. This index reflects fire risk and associated smoke-related air pollution, both of which are known to exacerbate respiratory conditions and potentially increase medication demand, particularly during hot and dry periods. The FWI was aggregated at the daily and weekly levels, and spatially usingthe mean. 3 Methods
Consideringthe complex interdependent relationships among various climate factors and their subsequent impact on the demand for respiratory medications, an exploratory analysis was necessary to identify relevant variables and appropriate lag structures. Based on the insights gained from this preliminary investigation, a set of predictive models incorporatingthe selected features and lag configurations wasthen employed to forecast drug demand. 3.1 Frequency-domain causality analysis and feature selection
The exploratory analysis began with unconditional and conditional Granger causality spectrato detect temporal dependenciesthat might be missed by time-domain methods, especially in cyclical data. This approach was motivated by periodic patterns observed in the Auto-Correlation Function (ACF) of drug demand and key environmental variables such astemperature and the FWI, highlightingthe relevance of frequency-domain analysis. The bootstrap procedure proposed by Farn`e and Montanari was used to assess whether observed relationships across frequencies were statistically significant. This approach is particularly advantageous compared to the parametric alternative proposed by Breitung and Candelon because it exhibits increasing power, approaching one, asthe process becomes non-stationary. This makes it especially appropriate given the structural changes observed in the time series under investigation. In the framework introduced by Farn`e and Montanari, the unconditional Granger causality spectrum quantifies how past values of one series help predict another across different frequencies. Specifically, for a bivariate stationary process ( X t , Y t ) [T] of length T, the unconditional Granger causality spectrum of the effect variable Y t with respect to the cause variable X t, at frequency ω, is defined as:
h Y Y ( ω ) h X→Y ( ω ) = log � σ 2 |P [˜] Y Y ( ω ) |
, �
where h Y Y ( ω ) denotesthe spectral density of Y t at frequency ω, P [˜] Y Y ( ω ) isthe element of the transformed transfer function matrix obtained from the normalized VAR representation of the process, and σ 2 isthe variance of the innovation in Y t after applyingthe linear transformation that orthogonalizesthe innovations. A positive value of h X→Y ( ω ) impliesthat past values of X t improve the linear predictability of Y t at frequency ω, withthe corresponding cycle length given by ω [−] . For each of the N stationary bootstrap samples simulated from the observed time series ( X t , Y t ) following Politis and Romano, resampled versions X t [∗] [and] [ Y] t [ ∗] [are generated. A VAR model] is estimated on each bootstrap sample using Seemingly Unrelated Regressions (SURE).
The optimal model order is determined by minimizingthe BIC, which balances model parsimony and goodness of fit. Then, the Granger causality measure h Y ∗ →X ∗ at Fourier frequencies f i = iT [−], for i = 1 , . . ., ⌊T/ 2 ⌋, is computed for each bootstrap sample. The median of these causality measures across all bootstrap samples is used to find the quantile of the bootstrap median distribution. At each frequency f i, the observed causality is flagged as significant if it exceedsthe corresponding quantile of the bootstrap distribution. A conservative Bonferroni correction is applied to control for multiple testing acrossthe entire frequency domain. In this framework, the significance level for each individual frequency is adjusted to 2 αF [−], where F representsthe total number of frequencies under consideration. This adjustment controlsthe family-wise error rate, thereby ensuring valid inference when testing for Granger causality at multiple frequencies simultaneously. While the unconditional approach reveals overall dependencies, the conditional Granger causality Spectrum extendsthe unconditional version by controlling for an additional influencing variable, enablingthe isolation of the direct causal effect between drug demand and climate-related variables while accounting for a potential confounder Wt:
˜ h ˜ Y Y ( ω ) h X→Y |W ( ω ) = log � σ Y Y |Q Y Y ( ω ) |
, �
where Q Y Y ( ω ) isthe element of the transfer function matrix derived from the VAR representation of the residual processes obtained by projecting both X t and Y t onto Wt and removingthe fitted components. σ Y Y denotesthe variance of the residuals from the second equation of the conditional VAR, where Y t isthe dependent variable after controlling for Wt . The bootstrap methodology in this case requires estimating a VAR model on Y t and Wt, allowing for the simulation of ( Y t [∗] [, W] t [ ∗] [) via residual bootstrap. Subsequently, the series] [ X] t [∗] [is simulated independently] using a stationary bootstrap. Finally, the conditional Granger causality spectrum is computed on each simulated triplet followingthe same steps as in the unconditional procedure. Building on the insights gained from the frequency-domain analysis, the most influential variables for forecasting demand for respiratory medications were identified through feature importance analysis of the Random Forest model with Moving Block Bootstrap (MBB-RF). In standard Random Forests, each decision tree istrained on a bootstrap sample drawn with replacement from the original dataset, enhancing robustness by reducing overfitting. However, this approach assumes independent observations, making it unsuitable for time series data, whichtypically exhibit temporal dependence. This limitation can be overcome throughthe use of the Moving Block Bootstrap (MBB) approach, as defined by K¨unsch, which samples consecutive blocks of data, preservingthe temporal structure and ensuringthat the inherent time dependencies are maintained. The resulting approach, here referred to as Moving Block Bootstrap Random Forest (MBB-RF), was previously considered by Goehry et al. as an adaptation of the Random Forest methodology to time series data. The feature importance metric adopted for feature selection isthe impurity-based measure, which quantifiesthe average reduction in residual sum of squares across all splits in the forest where the given predictor is used. In parallel to the nonparametric analysis performed throughthe MBB-RF, a sparse VAR model was also estimated to gain a complementary understanding of the dynamic interrelationships between drug demand and environmental predictors. Specifically, we employed a penalized VAR model with LASSO (Least Absolute Shrinkage and Selection Operator) penalty. This approach simultaneously enables variable selection and parameter estimation, making it well suited for high-dimensional time series with potentially correlated regressors. The use of a sparse VAR provides an interpretable framework in which only the most relevant temporal dependencies are retained, with many coefficients shrunk exactly to zero.
Algorithm 1 MMB-RF: Random Forest Regression with Moving Block Bootstrap Aggregation Input: Time series data D = { ( X t , y t ) } [T] t =1 [, where] [ X] [t] [ denotesthe vector of predictors at time] [ t] [, including] lagged values of boththe target variable and the covariates. B : Number of trees in the ensemble m : Number of variables randomly selected at each split n min : Minimum node size (stopping criterion)
1: for b = 1 to B do
2: Bootstrap: Generate a resampled dataset D [∗] [(] [b] [)] as follows:
(a) Define T − ℓ + 1 overlapping blocks of length ℓ :
Bi = { ( Xi , y i ) , .
. ., }, i = 1 , . . ., T − ℓ + 1
(b) Draw K = ⌈T/ℓ⌉ blocks with replacement from {Bi } (c) Concatenate the sampled blocks in temporal order:
D [∗] [(] [b] [)] = concat
(d) Truncate D [∗] [(] [b] [)] to length T if necessary
3: Tree Growing: Train regression tree Tb on D [∗] [(] [b] [)]
4: while node size > n min and variance > 0 do
5: Randomly select m predictors
6: Split the node usingthe feature and threshold that minimize the within-node variance
7: Grow child nodes recursively
8: end while
9: end for
Prediction: For anew observation x,
y ˆ RF ( x ) =
B
B � Tb ( x )
b =1
3.2 Forecasting models
Establishing a solid forecasting baseline began with implementingthe Prophet model, which was selected for its ability to capture complex temporal dynamics, including multiple seasonalities and structural breaks. This choice is particularly appropriate given that the time series of prescription pharmaceutical sales shows seasonal behavior and structural changes in the level.
The Prophet model, which relies solely on the target’s past values, offering a purely autoregressive benchmark of the target series, serves as a benchmark to assess whether the inclusion of climate covariates significantly improves forecast accuracy. Let S be the number of changepoints and at a binary vector indicating whether time t occurs after each changepoint s j, for j = 1 , . . ., S . Prophet decomposesthe observed time series y t into an additive model includingtrend ( gt ), seasonality ( st ), holidays ( ht ), and a Gaussian error term ϵ t ∼N:
y t = gt + st + ht + ϵ t . The trend component gt is modeled as a piecewise linear function:
gt = k + a [⊤] t [δ] t + m + a [⊤] t [γ] , � � � �
where k isthe base growth rate and δ adjuststhe rate at each changepoint. The offset term includesthe initial offset m and a correction vector γ, where each γ j = −s j δ j ensures continuity at changepoints. Prophet automatically selects changepoints from a predefined grid by placing a sparse prior on δ . The seasonal component st is expressed using atruncated Fourier seriesto provide a flexible model of periodic effects:
2 πnt � + b n sin � P
, ��
st =
N �
n =1
2 πnt � an cos � P
where P denotesthe period of seasonality and N controlsthe number of harmonics. Finally, the holiday component ht capturesthe effect of known recurring events by introducing indicator variables for each holiday and estimatingtheir impact with separate parameters. A VARX model, a MBB-RF, and a LSTM neural network were employed to incorporate climaterelated variables into the prediction of drug demand. These models were selected to represent a diverse set of methodological approaches—linear, ensemble-based, and deep learning—each capable of capturing different aspects of the relationship between climatic conditions and respiratory drug consumption. After verifyingthe absence of cointegration amongthe selected variables usingthe Johansen test, a VARX model with exogenous regressors was estimated. Considering a set of K endogenous variables and M exogenous predictors, a Vector Autoregressive model with exogenous regressors (VARX) can be expressed as y t = ν +
p � A i y t−i + Bx t + u t ,
i =1
where y t is a K -dimensional vector of endogenous variables, x t is an M -dimensional vector of time-aligned exogenous covariates, and ν is a K -dimensional vector ofintercepts. The matrices A i ( K × K ) capture the autoregressive structure over p lags, while B ( K × M ) quantifiesthe contemporaneous effects of the exogenous covariates. The innovation term u t is assumed to be a white noise process with zero mean and constant covariance matrix. The optimal lag order pis selected usingthe BIC. Althoughthe marginal distributions of the endogenous variables exhibit limited kurtosis, several display noticeable skewness or bimodality. Furthermore, Shapiro–Wilk tests reject the null hypothesis of normality for most variables at the conventional 5% significance level. The VARX model was estimated using a residual bootstrap procedure to ensure valid inference in the presence of such deviations from normality. This approach does not rely on asymptotic distributional assumptions and yields more accurate confidence intervals and test statistics in finite samples.
Furthermore, it enhancesthe robustness ofinference even when residuals appear approximately normal, as it accommodates nonstandard features in the marginal distributions of the input variables. Centered residuals ˆ u t = u t − u ¯ from are sampled with replacement to obtain bootstrap residuals u [∗] t [, which are then used to generate bootstrap samples] [ y] t [∗] [(] [b] [)] recursively, for b = 1 , . . ., B :
y t [∗] [(] [b] [)] = ˆ ν +
p � A ˆ i y t [∗] − [(] [b] i [)] [+ ˆ] [Bx] [t] [ +] [ u] [∗] t [(] [b] [)] , t = p + 1 , . . ., T. i =1
For each simulated series y t [∗] [(] [b] [)], bootstrap estimates A [ˆ] [(] i [b] [)] and B [ˆ] [(] [b] [)] are computed. Bias-corrected coefficients are then obtained as
�
�
A ˜ i = ˆ A i −
1
B �
B
ˆ
� A [(] i [b] [)] − A [ˆ] i b =1
, B ˜ = ˆ B −
1
B �
B
ˆ ˆ
� B [(] [b] [)] − B
b =1
. The MBB-RF approach was employed as a machine learning forecasting method due to its ability to capture complex nonlinear relationships and interactions among predictors without requiring strong parametric assumptions. The training procedure adopted for forecasting mirrorsthe one used duringthe variable selection presented in the previous subsection, ensuring methodological consistency in the modeling of temporal dependencies and lagged structures. A Long Short-Term Memory (LSTM) neural network was implemented to assess whether a deep learning model offers performance gains over traditional approaches in capturingthe potentially nonlinear and complex relationship between climatic factors and drug demand.
LSTMs are specifically designed for sequential data and are capable of learning long-range dependenciesthrough memory cellsthat maintain internal states over time. At eachtime step t = 1 , . . ., T, the state s c j ( t ) of memory cell c j evolvesthroughthe interaction ofinput and output gates, which regulate the flow ofinformation. This gating mechanism allowsthe network to retain or discard information selectively, thus effectively mitigatingthe vanishing gradient problem that affects standard Recurrent Neural Networks (RNNs) and limitingtheir ability to model longterm dependencies. The update rule for the internal state is:
�
s c j ( t ) = s c j + y [in] [j] ( t ) · g
�� u w c j ,u y [u]
u
, s c j = 0 ,
where g ( · ) is anonlinear activation function that controls how newinformation is incorporated into the memory and y [u] denotesthe output of unit u at the previoustime step. The input and output gate activations are defined as:
�
�
y [in] [j] ( t ) = fin j
�� u win j ,u · y [u]
u
�� u w out j ,u · y [u]
u
; y [out] [j] ( t ) = f out j
,
where win j ,u , w out j ,u ∈ Rare trainable connection weights determiningthe influence ofinput u on the respective gate unit j and fin j and f out j are activation functionsthat control the openness of the input and output gates. The output of the memory cell isthen given by:
y [c] [j] ( t ) = y [out] [j] ( t ) · h ( s c j ( t )) ,
where h ( · ) is an output activation function applied to the internal state. 5 If the corrected parameters ˜ A i and ˜ B violate the stability condition (i.e., at least one eigenvalue of the companion matrix lies within or on the unit circle), the biasterms are progressively shrunk by a factor δ < 1 until stability is restored. SHapley Additive exPlanations (SHAP) values were employed to enhance the LSTM model interpretability and elucidate the contribution ofindividual predictors. SHAP is grounded in cooperative game theory and assigns each feature an importance score correspondingto its average marginal contribution across all possible subsets of features. This yields explanationsthat are both locally accurate and globally consistent, effectively integratingthe strengths of global attribution methods, such as Permutation Importance, and local interpretability tools like LIME. Differently from Permutation Importance, SHAP is robust to multicollinearity, as it evaluates conditional expectations rather than marginal perturbations. Moreover, it avoidsthe instability and inconsistency that may affect LIME’s local surrogate models. Despite its computational cost, SHAP’s additive nature and sound theoretical grounding make it especially well-suited for interpreting complex, nonlinear architectures such as LSTMs. 3.3 Performance evaluation criteria
Forecasting models were trained on the first 338 weekly observations (approximately six and a half years) and evaluated on a hold-out test set consisting of the final 52 weeks.
A set of widely used error metrics was computed on the test set to evaluate forecasting performance, includingthe Mean Absolute Percentage Error (MAPE), Root Mean Squared Error (RMSE), RMSE-observations Standard deviation Ratio (RSR), and the seasonal Mean Absolute Scaled Error (MASE). Let y h denote the observed value at horizon step h, and ˆ y hthe corresponding predicted value, for h = 1 , . . ., H, where H = 52 isthe test set length. The error metrics are formally defined as:
~~�~~ ~~�~~ �
H
�
H
ˆ
�( y h − y h ) ,
h =1
 , RMSE = �
MAPE =
H
H �
h =1
ˆ y h − y h � y h
H �
1
H
H
ˆ
� |y h − y h |
h =1
RSR =
~~�~~ ~~�~~ � � � � � � � �
H , MASE = �( y h − y ¯ train )
h =1
H �
ˆ
�( y h − y h )
h =1
. T � |y h − y t−m |
h = m +1
1
T −m
T �
Here, T isthe length of the training set, ¯ y train denotesthe mean of the observed values in the training set, and m representsthe seasonal lag used for scaling in the MASE metric, which is set to m = 52 to reflect the annual seasonality inherent in weekly data. The use of multiple evaluation metrics enables a more nuanced understanding of model performance, capturing scale-independent errors, sensitivity to outliers, and relevance to domain-specific seasonal patterns. Specifically, MAPE assigns equal percentage weight to each observation, making it useful for evaluating relative errors across scales. RMSE quantifiesthe magnitude of typical prediction errors, placing greater emphasis on large deviations due to the squaring of residuals. Notably, RSR is mathematically related to R viathe identity R = 1 − RSR, allowing for a complementary interpretation. RSR provides anormalized error version of RMSE, facilitating comparisons across series with different levels of variability, while R measuresthe proportion of variance in the observed series explained by the model, making it particularly useful for assessing overall model fit. Finally, the seasonal version of MASE evaluates forecast accuracy relative to ana¨ıve yearly benchmark, making it particularly suitable for time series with recurring seasonal
patterns. Each metric was calculated individually for the forecasting models, allowing for a detailed assessment of predictive accuracy and a direct comparison of their effectiveness in forecasting drug demand. This evaluation also offered insight into the added value ofincorporating climate-related variables by benchmarking each model against a univariate Prophet baseline relying solely on the autoregressive dynamics of the target series. It is important to note that amongthe implemented models, the VARX model is distinct in its forecasting strategy, as it operatesthrough a recursive one-step-ahead approach, in which each prediction is iteratively fed back into the model to generate the subsequent forecast. In contrast, the Prophet, MBB-RF, and LSTM models are trained to produce the entire 52-week forecast horizon in a single step. This distinction is particularly relevant as recursive approaches in autoregressive frameworks are prone to error accumulation, particularly over longer horizons. Conversely, direct forecasting strategies mitigate the propagation of errors but may require larger training datasetsto achieve reliable performance. This comprehensive evaluation framework ensuresthat model comparisons are both methodologically rigorous and practically informative, particularly in guiding policy decisions related to respiratory health planning. 3.4 Software and implementation
Forecasting models such as Prophet, VAR, and MBB-RF were implemented in `R` usingthe `prophet`, `vars`, and `rangerts` packages, respectively. The LSTM neural network was developed in `Python` withthe `TensorFlow` and `Keras` libraries. Data preprocessing, visualization, and performance evaluation were primarily conducted in `R`, complemented by `pandas`, `matplotlib`, and `numpy` for auxiliary tasks in Python. Frequency-domain Granger causality tests, includingthe conditional version, were carried out in `R` viathe `grangers` package, which implementsthe methodology of Farn`e and Montanari. Additionally, the `sparsevar` package was employed to explore sparse VAR representations and dynamic interactions in high-dimensional settings. Interpretability of the LSTM model was addressed in `Python` usingthe `shap` library to compute Shapley values and evaluate feature contributions. Predictive uncertainty was assessed through Monte Carlo dropout by generating multiple stochastic forward passes during inference.
Meanwhile, the spectra for extreme rainfall and temperature standard deviation are comparatively weak, with significant peaks restricted to anarrow range of frequencies—high for extreme rainfall and low for temperature variability. These patterns suggest that, although extreme rainfall and temperature variability may have some localized, short-lived impacts, they are not consistent or dominant drivers of respiratory drug demand in the longer term. This highlightsthe importance of distinguishing between transient weather anomalies and systematic climatic shifts, especially in the context of climate change, which is expected to increase boththe frequency and intensity of extreme events. Conditional Granger causality spectra offer deeper insights into the distinct role of each climate variable by controlling for the confounding influence of the others, enablingthe disentanglement of direct from indirect effects and assessthe unique contribution of each variable to the prediction of drug demand. Consistent withthe unconditional results, temperature, specific humidity, wind speed, and the FWI remain amongthe most robust predictors, thoughthe signal for wind speed appears comparatively weaker. Their conditional spectra frequently exhibit significant peaks, especially at low frequencies, even after adjusting for other variables. This pattern suggeststhat part of their predictive power is shared with other climatic factors, yet their influence is not fully mediated, reinforcingtheir direct relevance for respiratory drug demand. Certain interactions between variables reveal more intricate dynamics. Notably, using precipitation as a conditioning variable enhancesthe spectral intensity of otherwise moderate predictors such as wet days, cloud cover, and extreme rainfall. Likewise, the spectrum of precipitation itself becomes more pronounced when conditioned on wet days or extreme rainfall. These mutual amplifications likely reflect shared meteorological processes—such asthe temporal clustering and persistence of precipitation events—that emerge more clearly once redundant variability is partialled out. By contrast, specific humidity as a conditioning variable tendsto flatten most spectra, substantially reducingthe apparent influence of many predictors. Even the spectrum of temperature, typically a dominant climatic driver, becomes less pronounced. While still statistically significant, it shows a flatter shape, with a subtle and low peak at low frequencies. This likely reflectsthe strong correlation between humidity and temperature, particularly in their seasonal behavior, suggestingthat humidity absorbs a considerable portion of the shared explanatory power when both are included. The behavior of FWI is also informative. Its conditional spectra remain strong across most combinations, reaffirmingthe robustness observed in the unconditional analysis. However, its predictive signal weakens considerably when temperature or humidity are used as conditioning 0.3
0.2
0.1
0.0
0.06
0.04
0.02
0.00
0.04
0.03
0.02
0.01
Unconditional Granger Causality Spectra: Climate Variablesto Drug Demand
Cloud Cover Extreme Rainfalls FWI
0.03
0.3
0.02
0.2
0.01
0.1
0.0
Precipitation Relative Humidity Temperature
2.0
1.0
1.5
1.0
0.5
0.5
0.0
0.0
Temperature SD Wet Days Wind Speed
0.08
0.09
0.06
0.06
0.04
0.02
0.03
0.00
0.00
0.0 0.1 0.2 0.3 0.4 0.5 0.0 0.1 0.2 0.3 0.4 0.5 0.0 0.1 0.2 0.3 0.4 0.5
Fourier frequency
Figure 2: Unconditional Granger causality spectra from climatic variablesto drug demand. Each panel displaysthe frequency-specific Granger causality spectrum from a given predictor to the weekly quantity of respiratory drug prescriptions. The solid line representsthe estimated spectrum. The dashed line denotesthe 95% significance threshold under the null hypothesis of no Granger causality, based on the empirical distribution of the bootstrap test statistic.
The dotted line correspondsto the Bonferroni-adjusted threshold that controls for multiple testing across frequencies. variables, indicatingthat a substantial part of FWI predictive power may be mediated throughthese two core drivers of atmospheric and fire-related dynamics. This interpretation aligns withthe known formulation of FWI, which includes components such as air temperature, specific humidity, and wind speed. Interestingly, temperature and FWI as conditioning variables enhance the spectral signatures of several weaker predictors, including wet days, cloud cover, and temperature standard deviation.
These patterns suggest that once the dominant seasonal component of temperature or fireweather-related variability is controlled for, subtler dynamics in the remaining variables emerge. Nonetheless, temperature and humidity as causal variables lose spectral strength when FWI isthe ExtremeRainfalls FWI Precipitation SpecificHumidity Temperature
0.20
0.15
0.10
0.05
0.00
0.10
0.05
0.00
0.75
0.50
0.25
0.00
0.125
0.100
0.075
0.050
0.025
0.000
1.0
0.5
0.0
2
1
0
0.12
0.09
0.06
0.03
0.00
0.15
0.10
0.05
0.00
0.12
0.09
0.06
0.03
0.00
0.0 0.1 0.2 0.3 0.4 0.5 0.0 0.1 0.2 0.3 0.4 0.5 0.0 0.1 0.2 0.3 0.4 0.5 0.0 0.1 0.2 0.3 0.4 0.5 0.0 0.1 0.2 0.3 0.4 0.5
Fourier frequency
Figure 3: Conditional Granger causality spectra. Each panel showsthe frequency-specific conditional Granger causality from a climatic variable X t in row to drug demand, while controlling for athird variable Wt in column. The solid line representsthe estimated spectrum. The dashed line indicatesthe 95% significance threshold derived from bootstrap inference under the null hypoth esis of no conditional Granger causality. The dotted line showsthe Bonferroni-adjusted threshold, accounting for multiple testing across frequencies. conditioning variable, again indicating a high degree of shared variance likely tied to seasonal cycles. Overall, these findings underscore the importance of addressing multicollinearity and indirect causal pathways in climate–health models. They also highlight that even predictors with weak unconditional effects may become relevant under specific conditional scenarios, reinforcingthe utility of frequency-domain Granger causality in revealing hidden structures in environmental influences on health outcomes. CloudCover TemperatureSD WetDays WindSpeed
0.3
0.2
0.1
0.0
0.15
0.10
0.05
0.00
0.9
0.6
0.3
0.0
0.20
0.15
0.10
0.05
0.00
1.0
0.5
0.0
2.0
1.5
1.0
0.5
0.0
0.15
0.10
0.05
0.00
0.20
0.15
0.10
0.05
0.00
0.08
0.06
0.04
0.02
0.00
0.0 0.1 0.2 0.3 0.4 0.5 0.0 0.1 0.2 0.3 0.4 0.5 0.0 0.1 0.2 0.3 0.4 0.5 0.0 0.1 0.2 0.3 0.4 0.5
Fourier frequency
Figure 4: Conditional Granger causality spectra. Each panel showsthe frequency-specific conditional Granger causality from a climatic variable X t in row to drug demand, while controlling for athird variable Wt in column. The solid line representsthe estimated spectrum. The dashed line indicatesthe 95% significance threshold derived from bootstrap inference under the null hypoth esis of no conditional Granger causality.
The dotted line showsthe Bonferroni-adjusted threshold, accounting for multiple testing across frequencies. 4.2 Variable selection
The role of climate-related variables in predicting drug demand is further assessed through a variable importance analysis conducted usingthe MBB-RF model, as described in Algorithm 1. Since random forests were not originally designed for time series applications, it is necessary to define the appropriate number of lagsto include in the feature set prior to model training. The optimal lag order was identified usingthe BIC applied to VAR models estimated for both unconditional and conditional Granger causality testing. Since no model selected an order greater than four, this lag structure was uniformly adopted across all variables duringthe training phase. The weekly resolution of the data further supportsthis choice, as a lag of four weeks approximately correspondsto one month, aligning with common pharmaceutical purchasing behaviors. Indeed, medications prescribed for respiratory conditions are typically issued monthly due to clinical factors—such as standard dosing regimens and regular follow-up requirements—and administrative constraints, including limits on repeat prescription durations and the implementation of monthly therapeutic plans. Also, settingthe number of lagsto four satisfies a principle of model parsimony, limiting complexity without sacrificing interpretability. Each series was assessed for weak stationarity to validate the MBB approach. The RF model wastrained using 1,000 trees, each built on an MBB sample with a block size of 52 weeks, reflectingthe seasonal nature of the data. Since respiratory drug consumption follows a clear annual cycle, this block size preservesthe intra-annual temporal dependencies critical to model performance. Eachtree was grown using a minimum node size of five, following standard recommendations in machine learning literature for regression tasks. This setting balances model flexibility and overfitting control, especially in time series contexts where deeper trees may lead to instability due to autocorrelated inputs. Given the length of the training set, this setup yields 287 overlapping blocks of 52 weeks each. For eachtree, 7 blocks were drawn with replacement, resulting in atraining sample of 364 observations per tree. By concatenatingthe selected blocks, the resampled series retains itstemporal coherence within each block, ensuringthat key seasonal patterns are not disrupted duringthe bootstrap resampling process. The results from the Granger causality analysis were integrated to the variable importance rankings obtained from the MBB-RF procedure to guide the selection of predictors for the forecasting models. Given the strong autoregressive nature of the drug demand series, lagged values of the target variable were included as predictors in all forecasting models. Indeed, autoregressive terms emerge amongthe most important features, highlightingthe persistence and temporal structure of respiratory drug consumption. The resulting importance plot reveals a pronounced elbow shape, suggestingthat only a small subset of variables substantially contributesto predictive accuracy. This curvature is used to define athreshold for variable selection, effectively identifyingthe most informative predictors while discarding less relevant ones. The lagged drug demand variables dominate the top ranks along withtemperature lags. In contrast, extreme rainfall and temperature standard deviation appear amongthe least relevant predictors, which is in line withtheir limited significance in both unconditional and conditional Granger causality spectra. Amongthe climate variables, temperature and specific humidity consistently rank highly. However, their strong static correlation indicates substantial redundancy and raises concerns about multicollinearity, which is particularly problematic for VAR models. While MBB-RF and LSTM tend to be more robust to collinearity, the same predictor set was initially applied across all modelsto ensure comparability of forecasting performance. Despite the relevance of specific humidity, supported by both Granger causality and feature importance analyses, its inclusion alongside temperature introduced overlap without leadingto consistent improvements in accuracy. Therefore, specific humidity was excluded in favor of temperature, which empirically emerged asthe more informative and impactful predictor. This selection strategy, guided by the elbowin the importance curve, led to a parsimonious yet effective predictor set, composed solely of temperature and the lagged values of respiratory drug demand, capturing boththe intrinsic temporal dynamics of the target series and the influence of the most relevant climate-related driver.
A sparse VAR model with a LASSO penalty was estimated to further assessthe system’s predictive structure. Consistent withthe MBB-RF model, the lag order was DrugDemand Temperature Temperature Temperature
DrugDemand Temperature
SpecificHumidity SpecificHumidity SpecificHumidity DrugDemand DrugDemand FWI FWI FWI
SpecificHumidity CloudCover FWI
CloudCover WindSpeed CloudCover CloudCover
Precipitation WindSpeed
TemperatureSD Precipitation WetDays
TemperatureSD WindSpeed WetDays
TemperatureSD Precipitation WindSpeed WetDays
TemperatureSD Precipitation WetDays
ExtremeRainfalls ExtremeRainfalls ExtremeRainfalls ExtremeRainfalls
MBB RF Feature Importance for variable selection
0 10,000,000,000 20,000,000,000 30,000,000,000 40,000,000,000 Variance reduction
Figure 5: Feature importance scores computed duringthe variable selection phase using a Random Forest model with Moving Block Bootstrap (MBB-RF). Importance is quantified astotal impurity reduction, averaged over bootstrap replicates. The plot highlightsthe most informative predictors for inclusion in the forecasting models. set to p = 4 to capture monthly cyclic dynamics, as suggested by the autocorrelation functions. All variables were standardized to ensure a fair regularization process, preventing dominance by high-variance predictors. The resulting model is highly interpretable due to the sparsity induced by LASSO. In the equation for drug demand, the response exhibits strong autoregressive behavior, with a dominant lag-1 coefficient and smaller positive effects at lags 2–4. This confirmsthe temporal persistence supported by boththe Granger spectral analysis and the MBB-RF model. Among climatic variables, temperature displays consistently negative coefficients at lags 1 and 3, suggestingthat warmer weather reduces respiratory drug demand—likely due to lower infection transmission and improved air quality. This interpretation aligns withthe importance scores from the MBB-RF model, where temperature emerged as a key predictor. Other variables, such as wind speed and precipitation, show smaller and more delayed effects. Cloud cover appearsto have limited influence. Specific humidity was retained in the estimation for completeness but resulted in zero coefficients due to its high correlation withtemperature, consistent withthe variable selection in the MBB-RF framework. Overall, boththe sparse VAR and MBB-RF models support a common narrative: drug demand is primarily shaped by its own past values and lagged temperature effects, with secondary contri butions from other environmental factors. The convergence of results across linear (sparse VAR) and nonlinear (MBB-RF) models strengthens confidence in the robustness and generalizability of the identified predictors. Table 1: Non-zero coefficients for the drug demand equation in sparse VAR of order 4. Variable Lag 1 Lag 2 Lag 3 Lag 4
Temperature -0.1841 0 -0.1167 0
FWI -0.0251 0 0 0
WindSpeed 0 0 0.0174 -0.0018 ExtremeRainfalls 0 0 0 0
WetDays 0 0 0 0 TemperatureSD 0 0 0 -0.0201 SpecificHumidity 0 0 0 0 DrugDemand 0.2575 0.0139 0.1282 0.1433 CloudCover 0 0.0257 0 0
Precipitation 0 0.0024 0 0
4.3 Forecasting models
The forecasts produced by the univariate Prophet model serve as a benchmark for evaluating predictive performance, while the fitted values—representing smoothed versions of the endogenous variables—are included as additional covariates in all subsequent forecasting models. This strategy enablesthe modelsto leverage Prophet capacity to flexibly extract trend and seasonal components without modifyingtheir underlying parametric structure. A rolling-window forecasting experiment further supported the inclusion of fitted series, confirmingtheir contribution to improved forecast accuracy. It is noteworthy that the Prophet model did not incorporate holiday effects as external regressors.
An in-depth analysis revealed no statistically significant difference in drug demand between holiday and non-holiday weeks, thus justifyingtheir exclusion. VARX MODEL
In the multivariate VARX model, where the lag order was selected based on the BIC and the estimation was carried out using a residual bootstrap procedure to ensure robust inference, the endogenous variables comprise temperature and respiratory drug demand. This specification is supported by the combined evidence from Granger causality spectra, variable importance analysis, and the sparse VAR model, all of which consistently identified temperature asthe most influential climate-related predictor, alongside the strong autoregressive dynamics of the target variable. A set of single-harmonic Fourier terms was included as exogenous regressorsto account for the dominant seasonal component observed in the weekly data. This allowed the model to capture annual periodicity explicitly without requiring a higher lag order, thereby reducing model complexity while preserving explanatory power. Additionally, a binary dummy variable was introduced to account for the systematic drop in pharmaceutical sales duringthe week of August 15th—anational holiday period characterized by reduced healthcare activity, hightemperatures, and consistently the lowest demand observed throughout the year. Fitted values from univariate Prophet models for both endogenous variables were also included as exogenous inputs. While the main goal isto forecast drug demand, incorporatingthe smoothed Prophet series for temperature helps improve its forecast within the VARX structure. This enhanced representation of temperature propagatesthroughthe system viathe multivariate dynamics of the model, ultimately yielding more accurate predictions of drug demand. The resulting hybrid strategy combines Prophet flexibility in capturingtrend and seasonal components withthe VARX ability to model lagged interdependencies, offering a robust and comprehensive framework for forecasting drug demand under climatic variability. The results in Table 2 indicate that higher temperatures in the four weeks preceding each observation are generally linked to a decrease in drug demand, although only the third lag demonstrates a statistically significant effect. This finding aligns withthe established seasonality of respiratory illnesses, whichtypically decline during warmer periods. Warmer weather is associated with improved air quality, increased outdoor activity, and lower transmission rates of respiratory pathogens. Conversely, colder temperatures can exacerbate respiratory conditions by increasingthe spread of viral infections, reducing mucociliary clearance, and triggering bronchoconstriction, particularly in vulnerable populations. These mechanisms help explain the heightened drug demand during colder months and support the observed negative association between temperature and respiratory drug consumption. However, using weekly data may obscure the short-term impacts of extreme temperature events, such as heat waves, which can have acute and nonlinear effects on vulnerable populations.
More detailed data (e.g., daily resolution) are therefore desirable to fully capture these dynamics and disentangle the effects of gradual seasonal changes from those of short-lived but intense climatic episodes, which are expected to become more frequent and severe in climate change. Accordingto the chosen estimation procedure, key assumptions regardingthe residuals of the VARX model were assessed using empirical p-value approaches based on bootstrapped statistics, ensuring robustness against finite-sample distortions and model misspecification. Specifically, residual diagnostics included a Portmanteau test for serial correlation and an ARCH-LM test for conditional heteroskedasticity. Bothtests were applied using 12 and 52 lags, corresponding approximately to a quarterly and an annual window, respectively, to capture potential short- to long-term volatility clustering in the weekly data.
Results were consistent across both lag specifications, showing no significant evidence of either serial correlation or heteroskedasticity, supportingthe validity of the residual bootstrap inference framework and the overall reliability of the estimated VARX model. Additionally, Granger causality tests in the time domain were conducted within the complete VARX system to evaluate the marginal contribution of each climatic predictor to drug demand over time. Althoughthe model structure differs from the spectral analysis due to the inclusion of exogenous regressors, the results provide strong evidence that temperature Granger-causes drug demand, with an empirical p-value of 0.0001 based on bootstrap inference. The p-value was computed using a conservative correction, which adds one to boththe numerator and denominator of the empirical ratio, thereby avoidingthe occurrence of zero p-values and yielding a more robust estimate in finite samples. While this confirms a highly significant predictive relationship, it does not provide any information on the frequency-specific structure of this link, which was previously revealed by the frequency-domain Granger causality spectra. The Impulse Response Functions (IRFs), displayed in Figure 6, offer valuable insights into the dynamic relationship between temperature, past drug demand, and current prescription patterns. These IRFs are derived from the estimated VARX model described in Table 2 to trace the effects of one-standard-deviation shocks on respiratory drug demand over time in each endogenous variable. The left panel in Figure 6 displaysthe response of drug demand to a one-standard-deviation shock in temperature. Althoughthe confidence bands encompass zero across much of the horizon—likely due to the inclusion of strong seasonal controls such as Fourier terms and the midAugust dummy—statistically significant responses emerge in the early weeks. Specifically, the effect initially turns negative, suggestingthat a sudden temperature increase may temporarily reduce demand, potentially due to behavioral adjustments (e.g., reduced indoor crowding or lower viral transmission in ventilated environments). This is followed, however, by a delayed increase in drug demand after the fourth week. The right panel in Figure 6 presentsthe IRF of drug demand to its own shock. The response confirmsthe strong autoregressive nature of the series, with an immediate and persistent effect that gradually diminishes over the 20-week horizon. This pattern is consistent with previous evidence from the MBB-RF importance analysis and aligns withthe structure of the VARX model, both of which highlight the predictive role of past demand levels. The IRF analysisthus addstemporal nuance to the interpretation of climate-health relationships, capturing both immediate and lagged effects within a multivariate framework. These results are further supported by the Forecast Error Variance Decomposition (FEVD) reported in Table S6 in the Supplementary material, quantifyingthe relative contribution of temperature shocksto the forecast error variance of drug demand. Althoughthe majority of forecast error variance is, as expected, attributed to the autoregressive dynamics of drug demand itself, temperature consistently explains around 6% of the variance from week 4 onward.
The decomposition stabilizes by week 12 and remains remarkably constant over longer horizons, with 95% confidence intervals for the temperature contribution ranging approximately from 2.2% to 10.4%. Random Forest
The same set of predictor variables was used across all modelsto ensure a fair comparison between the forecasting performance of the parametric VARX model and the nonparametric machine learning alternative represented by MBB-RF. Consistently withthe feature importance analysis, the predictive MBB-RF wastrained using 1,000 trees. Eachtree was built on a resampled time series composed of overlapping blocks of 52 consecutive weeks, thereby preservingthe seasonal structure inherent in the weekly data using a minimum node size of five. The resulting model achieved an out-of-bag (OOB) RMSE of 24,460.15, along with a RSR of 0.6033 and an R of 0.6360. The OOB metrics slightly outperform those obtained under a 52-week training–test split. This difference can be attributed to the nature of the MBB resampling scheme, which does not fully reproduce the temporal independence typically ensured by a strict separation between training and testing sets despite its statistical robustness. Specifically, the sampled blocks used to train eachtree remain close in time to the OOB observations, often sharing similar seasonal patterns. As a result, the forecastingtask becomes relatively easier than in a forward-chaining evaluation, where the model must extrapolate to unseen future periodsthat different environmental or behavioral dynamics may influence. In this context, test-set and rolling-window evaluation procedures offer a more conservative and realistic benchmark for assessing generalization performance under temporal dependence. Accordingto the feature importance results shown in Figure 7, the fitted values from Prophet rank at the top, suggestingthat the structural components extracted by this model provide substantial predictive value when used as covariates. These are followed by lagged values of drug demand and temperature, reinforcingthe findings from Granger causality analysis and autoregressive modeling. Notably, the first lag of temperature and the Prophet-fitted values outperform even the first lag of the target variable, despite the known Table 2: Bootstrap estimates and confidence intervals for the VARX model coefficients. Asterisks () indicate significance at the 5% level based on 1,000 bootstrap replications. Equation Coefficient Estimate Lower Bound Upper Bound
Drug Demand Constant -67722.6388 -153020 24610 Temperature -342.0322 -1939.9 1134.7 Drug Demand 0.1363 0.0135 0.2096 Temperature 887.7639 -803.26 2492.4 Drug Demand -0.0206 -0.1210 0.0650 Temperature -3013.8000 -4650.2 -1468.7 Drug Demand 0.1719 0.0561 0.2556 Temperature 920.8347 -653.29 2312.7 Drug Demand 0.2236 0.1199 0.3093 Prophet fitted (Drug Demand) 0.8742 0.6905 1.1375 August 15th -36094 -55318 -18324 Prophet fitted (Temp) 1106.3878 -2592.3 5188.7 sin -20919.3932 -40735 -653.71 cos -8007.0490 -43405 26181
Temperature Constant -4.3895 -9.9889 3.1394 Temperature 0.3308 0.1875 0.4078 Drug Demand 4.87e-6 -1.67e-6 1.24e-5 Temperature 0.0361 -0.1005 0.1302 Drug Demand -1.17e-6 -7.74e-6 5.92e-6 Temperature -0.0186 -0.1283 0.0852 Drug Demand -4.29e-6 -1.08e-5 2.09e-6 Temperature 0.0978 -0.0266 0.1755 Drug Demand 4.12e-6 -3.09e-6 1.05e-5 Prophet fitted (Drug Demand) 5.06e-7 -1.77e-5 1.61e-5 August 15th -0.1689 -1.6623 1.2397 Prophet fitted (Temp) 0.7814 0.5431 1.1008 sin 1.5482 -0.0743 2.5364 cos 1.4369 -1.1401 3.5149 IRF of Drug Demand on Drug Demand
0 10 20
Time
6000
3000
0
−3000
−6000
IRF of Temperature on Drug Demand
0 10 20
Time
25000
20000
15000
10000
5000
0
Figure 6: Impulse Response Functions (IRFs) estimated from a VARX model withtemperature and respiratory drug demand as endogenous variables. Exogenous regressors include a single harmonic Fourier term, Prophet-fitted values of the endogenous variables, and a dummy variable for the national holiday on August 15th.
This may reflect the ability of these inputsto capture broader, smoothed temporal signalsthat complement the more localized yet potentially noisier autoregressive dynamics. Seasonal regressors appear lower in the ranking but still retain non-negligible importance, confirmingthe role of annual cyclicality in shaping weekly drug demand. LSTM
The LSTM architecture was implemented on an input consisting of sequences of 52 time steps, each including drug demand and temperature. An extensive grid search over key architectural and training hyperparameters was performed to identify the optimal neural network configuration. Specifically, multiple combinations of the number of LSTM units, batch sizes, the inclusion of attention mechanisms, and early stopping strategies were considered.Each configuration wastrained usingthe Adam optimizer for up to 100 epochs, with early stopping based on the validation loss. A validation set comprising 10% of the training data was used to monitor performance, which was evaluated based on the minimum validation Mean Squared Error (MSE) achieved duringtraining. The best-performing architecture consisted of 80 units, a batch size 32, no attention mechanism, and early stopping enabled.
A 20% dropout rate was reactivated during inference using Monte Carlo Dropout to evaluate the stability and reliability of the LSTM forecasts. MBB RF Feature Importance for the predictive model
ProphetRxDemand
Temperature
DrugDemand
Temperature
Temperature
ProphetTemperature
DrugDemand
Temperature
DrugDemand
DrugDemand
CosFourier
SinFourier
August15th
0.0e+00 2.5e+10 5.0e+10 7.5e+10 1.0e+11 Importance
Figure 7: Feature importance derived from the predictive Random Forest model trained under a Moving Block Bootstrap scheme (MBB-RF).
Importance is measured by total impurity reduction across all trees. Variables are ranked in decreasing order of their contribution to forecasting weekly respiratory drug demand. The permutation importance analysis and SHAP values offer complementary insights into the relative influence of each predictor in the LSTM model. Across both interpretability frameworks, temperature emerges asthe most important variable. Its high permutation score indicatesthat randomly shufflingthis feature substantially degradesthe model’s accuracy, confirming its key role in shaping respiratory drug demand. The SHAP values corroborate this finding, showingthat temperature contributesthe largest marginal effect to model predictions acrossthe test set. The August 15th dummy and Prophet-fitted temperature are also relevant, particularly in the SHAP analysis. Their importance reflectsthe model sensitivity to known seasonal patterns and holiday-related drops in pharmaceutical activity. The Fourier terms exhibit moderate importance, suggestingthat while they explicitly encode seasonality, part of this information may already be captured by the LSTM internal state, which is inherently suited to learning recurringtemporal patterns from the input sequence. In contrast, the permutation analysis revealsthat Drug Demand and Prophet-fitted Drug Demand receive negative importance scores. This suggeststhat includingthese variables may introduce noise or redundancy, potentially degrading performance when not properly regularized. One explanation isthat the autoregressive signal in drug demand may already be sufficiently captured by the model temporal dynamics and temperature-driven effects, makingthe explicit inclusion of these features unnecessary or even detrimental. The SHAP values for these same predictors are relatively small, indicating a limited contribution to individual predictions and further supportingtheir marginal role in the model. Overall, the results suggest that environmental and seasonal exogenous features are more pivotal than the autoregressive component in shaping LSTM forecasts. This finding is particularly relevant in our interest in climate change, as it underscoresthe dominant influence of climatic conditions on drug demand. Nevertheless, to ensure a fair and consistent comparison across models, the same feature set was retained throughout all predictive experiments, thereby isolating model performance differences from input information differences. Discussion
Table 3: Forecasting performance of the four models—Prophet, VARX, MBB-RF, and LSTM—on one-year-ahead test data. Each model is estimated or trained on atraining set of 338 weekly observations and evaluated on a 52-week horizon. All models use the same predictors, including autoregressive lags, seasonal Fourier terms, the August 15th dummy, and Prophet-fitted temperature and drug demand values as exogenous regressors. Model MAPE RMSE RSR R MASE
Prophet 12.0000 36924.92 0.7395 0.4531 1.1207 VARX 10.8102 33243.92 0.6658 0.5567 1.0081
RF 10.5861 33361.86 0.6681 0.5536 0.9892
LSTM 11.3413 30886.91 0.6246 0.6099 0.8897
We assessthe one-year-ahead forecasting performance of four models using a harmonized set of predictors and a consistent rolling evaluation framework. All models were trained on a 338-week windowand evaluated over a 52-week forecast horizon. The input features include seasonal Fourier terms and Prophet-fitted values for temperature and drug demand.
This standardized setup ensuresthat differences in performance can be attributed to the model architecture itself, rather than to variation in input data or training conditions. Regarding relative accuracy, the MBB-RF model slightly outperforms all others, achievingthe lowest MAPE and a MASE below 1. These values indicate that MBBRF effectively captures average demand patterns and consistently improves upon a seasonal na¨ıve benchmark. Its R value confirms a solid proportion of variance explained, supporting its ability to generalize across diverse temporal contexts. The VARX model follows closely, with only marginally higher MAPE and MASE values. However, it achievesthe lowest RMSE and RSR, highlighting its effectiveness in reducing large forecast deviations—an important property when the cost of extreme prediction errors is high. Its R of 0.5567 isthe second-highest, reaffirmingthat VARX balances accuracy and model stability. Additionally, its structural interpretability offers analytical advantages, astechniques such as IRFs and FEVD allow for the examination of lagged dynamics and propagation mechanisms across variables. The LSTM model deliversthe best overall performance on several metrics: RMSE, MASE, and the highest R. These results underscore its capacity to learn longrange dependencies and model intricate temporal fluctuations. However, its MAPE is slightly higher than those of RF and VARX, suggesting atrade-off: while LSTM excels in minimizing squared error and explaining variance, it may be more sensitive to shorter training sequences or overfitting. Furthermore, its black-box nature limits interpretability, requiring post-hoc methods like SHAP or LIME to extract meaningful insights. Despite its proficiency in modelingtrends and seasonality, the Prophet model ranks lowest across all metrics.
With a MAPE of 12.00%, MASE of 1.1207, and a relatively low R, Prophet strugglesto matchthe performance of more flexible or feature-rich approaches. Its lack Figure 8: (a) Permutation importance
Figure 9: * (b) SHAP values
Figure 10: Model-agnostic interpretation of LSTM predictions. Panel (a) showsthe permutation importance scores, measured by the increase in MSE when each feature is randomly shuffled. Panel (b) displaysthe average SHAP values, quantifyingthe marginal contribution of each predictor to the model output acrossthe test set. of climatic covariates hampers responsivenessto environmental drivers, particularly critical in this application where temperature, the leading climatic driver, substantially improves performance, notably in RMSE, MASE, and R . Importantly, only MBB-RF and LSTM yield MASE values below 1, suggestingthat these models consistently outperform a seasonal na¨ıve benchmark. This reinforcestheir practical reliability, especially when evaluated using normalized accuracy metrics. In contrast, the slightly higher MASE of the VARX model—despite its strong RMSE and R scores—suggeststhat its forecasts do not consistently surpassthe accuracy of a simple seasonal repetition. This outcome may be attributed to a strongly seasonal demand pattern and limited abrupt climatic shifts within the observed period. In such cases, ana¨ıve seasonal forecast (e.g., using last year’s values) may perform reasonably well, setting a high bar for MASE-based comparisons. Nonetheless, VARX remains valuable for its ability to capture and interpret lagged dependencies and for offeringtransparent insights into longrun interactions, particularly when environmental variability is subtle and structural understanding is prioritized over marginal gains in relative accuracy. These findings underscore the importance ofintegrating climate-sensitive covariates into forecasting models. Approachesthat combine temporal dependencies with exogenous regressors, particularly those capable of modeling nonlinear interactions or complex sequences, consistently outperform simpler autoregressive alternatives. However, the robust performance of VARX demonstratesthat well-specified parametric structures can remain competitive, especially when model transparency and causal interpretability are essential. The final choice among models reflectstrade-offs among predictive accuracy, interpretability, and model complexity. This balance depends heavily on the domain context of forecast-driven decisions and policy relevance. When the primary objective is predictive accuracy, particularly in minimizing large forecast deviations, the LSTM model stands out asthe most effective choice. For applications where average predictive performance is prioritized, VARX and MBB-RF emerge as strong alternatives, achieving lower MAPE values. Amongthese, VARXis particularly well suited when interpretability and the analysis of temporal dynamics are essential, due to its parametric structure and the availability of tools such as impulse response functions and forecast error variance decomposition. In contrast, MBB-RF may be preferable in contexts where nonlinear relationships and complex feature interactions are expected, offering greater modeling flexibility. The final model selection should reflect the trade-offs between forecast precision, transparency, and the need for structural insight, depending on the specific goals of the study.
In summary, when accuracy under large deviations is paramount, LSTM is preferable; when average precision and transparency are valued, VARXis ideal; and in contexts where complex interactions may play a role, MBB-RF offers a flexible and robust alternative. 5 Concluding remarks
This study implements and compares a range of forecasting modelsto predict the weekly demand for prescription pharmaceuticals in Greece, withthe broader goal of quantifying how climate and environmental conditions influence medical needs in the population. As climate change intensifies—bringing more frequent and extreme weather events—anticipating shifts in drug demand becomes increasingly important for public health planning and supply chain resilience. We first apply a joint analysis of unconditional and conditional Granger causality spectrato understand which environmental factors matter most. This frequency-domain approach revealsthat temperature, specific humidity, wind speed, and the FWI exhibit the strongest and most consistent associations with respiratory pharmaceutical demand, particularly at lowand medium frequencies. Crucially, these variables maintain statistical significance even when controlling for other climate-related predictors, emerging as prime candidates for inclusion in forecasting models. However, the Granger causality framework inherently limitsthe analysisto pairwise or triplet relationships, making it difficult to assessthe joint influence of multiple predictors. To overcome this limitation, we complement the Granger-based findings with a Random Forest model trained using a Moving Block Bootstrap scheme (MBB-RF), which enables simultaneous evaluation of all predictors while accounting for temporal dependence. The MBB-RF importance rankings corroborate the Granger results: temperature, specific humidity, and the Fire Weather Index (FWI) consistently emerge asthe most influential predictors of pharmaceutical demand. The sparse VAR with Lasso penalization further supportsthe predominance of temperature, where the largest coefficients are associated with lags 1 and 3 of temperature, followed by lag 2 of cloud cover and lag 1 of FWI. This validates and reinforcesthe findings obtained throughthe other methodologies. Nevertheless, the importance distribution exhibits a pronounced elbow, indicating a sharp drop in variable relevance, measured in terms of impurity reduction, after accounting for temperature and the autoregressive terms of drug demand. Based on this observation, and to reduce redundancy, only temperature and the lagged values of the target variable were retained in the final predictor set. This parsimonious selection balances model simplicity and forecasting performance, capturing boththe intrinsic temporal dynamics of respiratory pharmaceutical use and the most impactful climatic signal. Performance metrics for the implemented forecasting models highlight that incorporatingtemperature as a climatic regressor consistently enhances forecasting accuracy across all model classes. Its seasonal behavior reflects not only climatic rhythms but also behavioral and epidemiological dynamics related to the spread of respiratory illnesses, making it an effective predictor of demand fluctuations. Amongthe evaluated models, the VARX structure that included Prophet-fitted exogenous inputs initially offered a good balance between interpretability and predictive accuracy based on squared error metrics. However, the MBB-RF outperformed it in terms of relative error measures such as MAPE and MASE, while also enabling variable importance analysisthat supported the findings from the Granger causality framework.
The LSTM model, although slightly less accurate in relative terms, achieved the best performance in RMSE and RSR, highlighting its strength in capturing nonlinear dynamics and high-frequency temporal patterns when sufficient structure is available in the data. Taken together, the evidence pointsto a clear and concerning conclusion: climatic conditions, especially temperature, have a causal impact on the demand for respiratory medication. As climate variability grows, this relationshipis likely to intensify, posing new challenges for healthcare systems. Forecasting models must therefore be equipped to anticipate these shifts, not only to optimize supply chains, but also to inform broader health preparedness strategies in a warming world. Another avenue for future development concernsthe temporal scope and resolution of the data. In this study, we relied on seven and a half years of weekly data, which is suitable for identifying recurring seasonal patterns and medium-term variability. However, to more robustly assessthe long-term effects of climate change, which unfold over decades, it would be valuable to extend the historical coverage of the dataset. Moreover, adopting a finer temporal resolution, such as daily data, could improve the model’s ability to capture the effects of short-lived but intense events, including heat waves, extreme rainfall, and other forms of acute environmental stressthat are becoming increasingly frequent due to climate change. This would enable a more detailed analysis of how climate extremes affect medical needs on shorter timescales, complementingthe broader patterns captured by weekly aggregates. Further research could also expand the analysisto include a spatial dimension, integrating geographic heterogeneity and localized climate exposure.
These enhancements could further support climate-aware planning in public health logistics, helping ensure that medication supply chains remain responsive in the face of a rapidly changing environment. Acknowledgments
The authors sincerely thank Dr. Filippo Benetti and Dr. Rosa Falotico of Alira Health Srl for their support in the acquisition and interpretation of the pharmaceutical data. The dataset used in this study is confidential and cannot be publicly shared. Author contributions
Viviana Schisa and Matteo Farn`e contributed to the conception and design of the work. Viviana Schisa contributed to the data management, statistical analysis and interpretation of results. Viviana Schisa and Matteo Farn`e contributed to the methodology and investigation. Viviana Schisa drafted the work, and all the authors reviewed it for important intellectual content. All the authors approved the publication of the final version of the manuscript.
All the authors agreed to be accountable for all aspects of the work to ensure that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. Financial disclosure
Viviana Schisa acknowledges financial support from Alira Healthand the European Union – NextGenerationEU, under the Italian National Recovery and Resilience Plan (PNRR), Mission 4, Component 2, Investment 3.3, as part of the research fellowship Ex D.M. 352/2022. Conflict ofinterest
The authors declare no potential conflict ofinterests. Supporting information
Additional materials supportingthe findings of this study are provided in the online Supplementary Information. These include detailed descriptive statistics and distributional diagnostics for all climatic and pharmaceutical variables used in the analysis, along with diagnostic plots such as density curves, boxplots, and autocorrelation functions. Univariate time series models are presented to characterize the seasonal and stochastic structure of respiratory drug demand, includingthe estimation of a SARIMA–eGARCH model with volatility and asymmetry diagnostics. Further sections provide extensive validation of the multivariate VARX framework, including normality and cointegration tests, rolling window forecasting results, and forecast accuracy comparisons across models. The impact ofincluding Prophet-fitted series as external regressors is evaluated, and additional analyses explore model robustness under temporal shifts. Finally, a map of data collection sites across Greece is included to document the spatial coverage of the dataset. Locations correspond to the most populous city in each administrative region, ensuring alignment with pharmaceutical consumption patterns.