Examining Forecasting Obstacles for Water Level in Lake Malawi

STATS 531/631 Midterm Report - Winter 2025

Introduction

Lake Chilwa, a 1,750 km² lake in Malawi, experiences unpredictable dry spells, occasionally drying out completely, located about 50 km east of Lake Malawi. This variability disrupts fisheries, threatens biodiversity, and reduces water availability for agriculture, posing challenges for local communities (Njaya et al. 2011). Furthermore, the difficulty in predicting its water levels makes it hard for residents to plan for these fluctuations, increasing their vulnerability to droughts and resource shortages (Chiwaula et al. 2012). Wishing to study this phenomenon but lacking sufficient data for Lake Chilwa, we turn our attention to Lake Malawi, hoping to uncover insights that may inform future research and management strategies for Lake Chilwa (Lancaster 1979).

Lake Malawi, Africa’s third-largest lake, is shared by Malawi, Mozambique, and Tanzania. It provides vital water for hydropower, farming, and household use (Mbaye et al. 2015). Covering 28,800 km² with a depth of up to 785 meters, the lake is home to many unique fish species, especially Cichlids (Mbaye et al. 2015). However, the lake’s water levels often change due to weather patterns and natural water movements. Rainfall, evaporation, and runoff are the main factors affecting the lake’s water levels (Kumambala 2010). Rain adds water, while evaporation, worsened by climate change, reduces it (Kumambala 2010). Runoff from the 126,500 km² surrounding area brings more water into the lake and the Shire River is the only outflow of the lake (Mtilatila et al. 2020). In the past, water levels have varied greatly, causing long droughts between 1915 and 1935 and severe floods in the late 1970s (Kumambala 2010),(Mtilatila et al. 2020).

This study will use data from the Climate Engine website¹ and the Malawi government² to explore which of the water level components: precipitation, evaporation, or runoff causes the greatest difficulty in predicting water levels. In water-scarce regions, simply confirming that water levels are difficult to predict may offer little practical value. However, because water levels function within a largely closed system, pinpointing which components contribute most to forecasting challenges can provide deeper insight into the underlying sources of uncertainty. This understanding may help improve water management strategies in areas facing water insecurity.

Subsequently, we seek to forecast these three variables and evaluate the quality of the forecast models produced. To contextualize the performance of these forecasts, it is important to consider their implications for predicting water level. By assessing both the predictive accuracy of a forecast and the relationship between the underlying variable and water level, we can better understand the extent to which forecast errors can impact water level projections.

Exploratory Data Analysis

Descriptive Statistics for Different Variables (mm)

Measure	Water Level	Evaporation	Precipitation	Natural Runoff
Mean	473472.99	85.44	98.03	27.85
Median	473485.00	75.59	36.51	2.22
SD	641.53	57.90	109.77	51.24
Minimum	472106.67	11.88	0.00	0.00
Maximum	475036.67	189.03	446.65	322.92
IQR	939.58	115.22	192.03	22.29
Skewness	0.02	0.23	0.81	2.34
Kurtosis	2.32	1.50	2.30	8.76

Lake Malawi Water Level Patterns (Graph 1)

The time series of Lake Malawi water levels from 1992 to 2017 reveals patterns of fluctuation driven by rainfall, evaporation, outflow, and natural runoff. Research suggests these fluctuations are significantly influenced by climatic variability, such as changes in rainfall patterns linked to phenomena like the El Niño-Southern Oscillation (Nicholson 2017). Peaks occur during periods of heavy rain or flooding, often associated with La Niña events, while declines occur with droughts, particularly during El Niño years. Additionally, high evaporation rates due to regional climatic conditions contribute to these trends (Kumambala 2010). The average water level of 473,473 mm, closely aligned with the median of 473,485 mm, suggests a symmetrical distribution with moderate variability (standard deviation of 641.53 mm). The interquartile range (IQR) of 939.58 mm highlights significant spread, while the low skewness (0.023) and moderate kurtosis (2.32) indicate a nearly normal distribution with occasional extremes.

Lake Malawi Precipitation Patterns (Graph 2)

Lake Malawi experiences distinct seasonal precipitation patterns, with rainfall peaking during wet months and dropping significantly in dry periods. From 1992 to 2017, the lake’s precipitation followed an annual cycle, influenced by broader climatic factors such as the Intertropical Convergence Zone (Nicholson 2017). The average monthly rainfall of 98.03 mm, with a median of 36.51 mm, suggests a right-skewed distribution where occasional heavy rainfall events contribute to overall precipitation. A standard deviation of 109.77 mm highlights the high variability, while a skewness of 0.81 and kurtosis of 2.30 indicate occasional extreme rainfall events (M. R. Jury and Mwafulirwa 2002). The interquartile range (IQR) of 192.03 mm further emphasizes the wide seasonal fluctuations in precipitation. This variability aligns with previous studies on Lake Malawi’s climate, which show that precipitation is a key driver of hydrological changes in the region (Bootsma and Hecky 2003).

Evaporation Patterns (Graph 3)

Lake Malawi’s evaporation follows a seasonal cycle, peaking in warmer months and declining in cooler periods due to temperature fluctuations. From 1992 to 2017, the mean evaporation of 85.44 mm closely aligns with the median of 75.59 mm, indicating a relatively balanced distribution. However, the standard deviation of 57.90 mm and an IQR of 115.22 mm highlight significant variability, with occasional high-evaporation months but few extreme outliers (skewness = 0.23, kurtosis = 1.50) (Bhave et al. 2020). Evaporation is a major contributor to water loss in Lake Malawi, often exceeding the Shire River outflow, making it a crucial factor in water balance modeling (Mtilatila et al. 2020). Climate change has intensified evaporation trends due to rising temperatures and changing wind patterns, further impacting water availability for hydropower, agriculture, and local communities (Lyons et al. 2015). These findings align with past studies on Malawi’s climate, which emphasize the growing importance of adaptive water management strategies to mitigate evaporation-driven declines in lake levels (Mark R. Jury 2014).

Natural Runoff Patterns (Graph 4)

The natural runoff time series from 1992 to 2017 shows repeating cycles, with peaks following heavy rainfall and low runoff during dry periods, reflecting regional rainfall patterns and influencing Lake Malawi’s water levels. The mean runoff of 27.85 mm is significantly higher than the median of 2.22 mm, indicating a highly skewed distribution where low runoff is common, but extreme runoff events occasionally occur (Bhave et al. 2020). The standard deviation of 51.24 mm and IQR of 22.29 mm highlight substantial variability, while a skewness of 2.34 and kurtosis of 8.76 confirm frequent extreme values. The minimum runoff is 0 mm, while the maximum reaches 322.92 mm, reflecting the sporadic nature of intense runoff periods (Mark R. Jury 2014). These extreme fluctuations directly impact Lake Malawi’s water balance, influencing both seasonal water levels and long-term hydrological trends. Increased runoff can lead to flooding, while prolonged low-runoff periods contribute to drought conditions, both of which pose challenges for hydropower generation and water resource management (Lyons et al. 2015).

Lake Malawi Water Level by Month (Graph 5)

Each month, Lake Malawi’s water levels fluctuate in response to precipitation, evaporation, and regulated outflows. During the rainy season (December–April), increased precipitation and runoff cause water levels to rise, peaking in March and April as the Intertropical Convergence Zone (ITCZ) shifts southward, bringing heavy rainfall (Nicholson 2017). As rainfall decreases in May, the lake enters the cool-dry season (May–August), where lower temperatures slow evaporation, but minimal inflows result in a gradual decline (Kumambala 2010). The most significant drop occurs between September and November, coinciding with peak evaporation rates, minimal precipitation (Bootsma and Jørgensen 2013). Climate variability further influences these patterns, with El Niño events reducing wet-season rainfall, leading to weaker peaks, while La Niña enhances precipitation, resulting in higher-than-average water levels (Mark R. Jury 2014). Rising temperatures are causing more evaporation, increasing water loss in the dry season and making the lake more vulnerable to climate changes (Lyons et al. 2015).

Model Fitting

In this section, we explore modeling the trajectory of evaporation, precipitation, and natural runoff in the Lake Malawi region. While we have data through August 2017, we will only use data through the end of 2013 for model building. We will then use data from 2014-2017 as test data to obtain error estimates for our forecasts, discussed more in the section below.

For illustrative purposes, we document the model-building for the evaporation data. We built models precipitation and natural runoff using similar model-building strategies.

Evaporation

We begin with a quick examination of the evaporation trend, shown in Graph 3 above. In this graph, we observe a clearly periodic relationship. As one may expect, we observe these annual oscillations for precipitation and runoff as well, with the periods corresponding to the seasonality of the water cycle (“wet/dry seasons”) (Chou et al. 2013). Subsequently, we consider seasonal ARMA models. As in class, we will begin by exploring models of the form \(SARMA(p,q)\times (1,0)_{12}\).³ We explore Akaike information criterion (AIC) values for such models in the table below (using `stats::arima()`` for model fitting).

stats::arima()-Produced AIC Table for Evaporation

	MA0	MA1	MA2	MA3	MA4
AR0	2151.35	2084.28	2072.52	2070.60	2071.23
AR1	2066.66	2068.36	2070.21	2071.67	2073.22
AR2	2068.32	2049.37	2051.01	2064.84	2059.49
AR3	2069.96	2051.04	2059.76	2062.26	2059.77
a All models contain a \((1,0)_{12}\) seasonal term.

In this table, we see some inconsistent AIC values. For example, the \(SARMA(2,2)\times(1,0)_{12}\) model fit has estimated AIC 2051.01, whereas the \(SARMA(2,3)\times(1,0)_{12}\) model fit has estimated AIC 2064.84. Assuming proper maximization, the latter should have an AIC of at most two more than the former. Given this is not the case, we conclude that stats::arima is not maximizing the likelihood properly. Subsequently, we turn to the arima2 package (Wheeler, McAllister, and Sylvertooth 2024).

arima2::arima()-Produced AIC Table for Evaporation

	MA0	MA1	MA2	MA3	MA4
AR0	2151.35	2084.28	2055.99	2070.60	2071.23
AR1	2066.66	2068.36	2070.21	2063.58	2058.96
AR2	2068.32	2049.37	2051.00	2007.06	2008.82
AR3	2069.96	2051.04	2045.77	2046.18	2010.31
a All models contain a \((1,0)_{12}\) seasonal term.

As before, we observe issues with AIC consistency. As a result, we turn our attention to the technical details of the likelihood optimization. We now fit SARMA models as follows:⁴

sarma_model_fit <- arima2::arima(x = time_series_data, order = c(p,0,q), 
                                 seasonal = list(order=c(P,0,Q), period=12),
                                 diffuseControl = FALSE, max_repeats = 25, 
                                 max_iters = 200, 
                                 optim.control = list(maxit=1000))

Here, setting diffuseControl=FALSE rejects the default method of ignoring the likelihood values of the initial observations if they’re controlled by the diffuse prior (i.e., have a Kalman gain of at least 1e4).⁵ We also use optim.control = list(maxit=1000) to allow each run of the optim() function to use more iterations to reach an estimated optimal value (default 100). We also increase the maximum number of random restarts for the CSS-ML method by setting max_iters=200 and we increase the threshold for the number of consecutive random restarts that don’t improve the likelihood by setting max_repeats=25.⁶

This generated the following table:

arima2::arima()-Produced AIC Table for Evaporation (With Adjusted Optimization Parameters)

	MA0	MA1	MA2	MA3	MA4
AR0	2151.35	2084.28	2072.52	2070.60	2071.23
AR1	2066.66	2068.36	2070.21	2063.58	2058.96
AR2	2068.32	2049.37	2051.00	2007.06	2008.82
AR3	2069.96	2051.04	2045.77	2046.18	2010.30
a All models contain a \((1,0)_{12}\) seasonal term.

Again, we observe some inconsistent AIC values - particularly surrounding the \(SARMA(2,3)\times(1,0)_{12}\) model. We investigate the (inverted) roots of this SARMA model below.

As we can see above, there are multiple roots on the unit circle (for both the AR and MA polynomials), likely contributing to the persistent numerical issues shown above. We explored simpler models (plots omitted for conciseness); however, we continued to observe a unit root for the MA polynomial whenever \(q>0\). Taking the lowest AIC of the models without MA unit roots resulted in the \(SARMA(1,0)\times(1,0)_{12}\) model (which did not have any roots on/in the unit circle). However, when we examine the sample autocorrelation of this model below, we observe a spike for 12 month lags.

This indicated that there was still seasonality in the residuals of the model. Therefore, we explored adding an MA term to the seasonal component of our SARMA model. I.e., we turned attention to models of the form \(SARMA(p,q)\times(1,1)_{12}\). Taking the model with the lowest AIC (among \(p\in\{0, ..., 3\}, q\in\{0,..., 4\})\), gave \(SARMA(1,2)\times(1,1)_{12}\).⁷⁸ We observe that this largely flattened the autocorrelation at lag 12. While there is some potential periodicity in the residual correlation, this has largely been reduced and we proceed.

While the fitted model appears to describe the data well, we observe that a Normal QQ plot shows highly non-normal behavior of the residuals of this model.

As suggested in class, we thus explore a log transformation. Thus, we consider models of the form \(SARMA(p,q)\times(1,1)_{12}\) for \(\log(Evaporation)\). Using AIC (with the diagnostic specifications discussed above) suggested that \(SARMA(2,0)\times(1,1)_{12}\) would be an appropriate model.⁹ We examine the sample ACF, fitted values and the normality of the residuals below.

As we can see above, the model fits the log evaporation data well and the sample autocorrelation shows little evidence of problematic autocorrelation in the residuals. Furthermore, while the QQ plot does not show perfect alignment with normal quantiles, we see that the log transformation has largely alleviated the non-normality in the residuals. Given these observations, we select this model for our time series analysis of evaporation.

Precipitation and Natural Runoff

We underwent a similar process as above for fitting models to the precipitation and runoff data. One key difference was that these time series had values of 0 for some months (i.e., \(precipitation_n, runoff_n\geq 0\), whereas \(evaporation_n>0\)). This prevented directly log-transforming the time series. As advised by (Hyndman), we explored two different transformations (both of the form \(\tilde{y}_n=\log(y_n+\lambda))\)):

• \(\lambda_1\): The square of the first quartile of the time series divided by the third quartile.
• \(\lambda_2\): One half the smallest non-zero value of the time series.

Having explored models under both transformations, we saw that the best-performing models (in line with the process we saw above for evaporation) for both precipitation and natural runoff used the \(\lambda_1\) transformation. This gave \(\lambda_{1,precip}\approx -0.00056, \lambda_{1,runoff}\approx 0.00003.\)

Using the same decision-making process as above, we modeled \(\log(Precip_n + \lambda_{1,precip})\) using a \(SARMA(3,1)\times(1,1)_{12}\) model, and \(\log(Runoff_n + \lambda_{1,runoff})\) using a \(SARMA(2,0)\times(1,1)_{12}\) model.

Model Evaluation

As discussed above, our primary research question involved analyzing whether any of evaporation, precipitation, or natural runoff were singularly difficult to predict (which would help explain the unpredictability in water level). We will use mean absolute error (MAE) as the primary metric to assess “unpredictability.” I.e., for each variable, we

• Calculate the residuals (in natural units)¹⁰ of the rolling one-step-ahead forecasts of the models above on the test data (i.e., measurements taken in 2014 and beyond).
• Take the average absolute residual over the test set.

Note that, by “rolling one-step-ahead forecasts,” we mean that we try to predict the next month’s evaporation/precipitation/runoff given all previous months’ measurements (although we do not re-fit any SARMA models in this process). Furthermore, by restricting this analysis to the test data (which was not used to train the model), the observed MAE should be a consistent estimator for \(\mathbb{E}\left[|\hat{Y}_n-Y_n|\right]\).

This analysis gives the following MAE values:

Mean Absolute Error (MAE) for Different Variables (mm)

Evaporation	Precipitation	Runoff
6.32	32.27	15.66

The table above suggests that evaporation is the most predictable of the three variables, whereas precipitation is the least. We see this in the plots below, as the fitted model has difficulties forecasting the magnitude of the seasonal peaks for precipitation.

Regression

As discussed above, we wish to contextualize our MAE estimates. For example, just because precipitation saw the highest MAE of the three variables considered, that does not immediately imply the inability to predict precipitation undermines the ability to predict water level. However, if precipitation were difficult to predict and precipitation was highly associated with water level, the evidence for such a claim strengthens.

In this section, we explore regression analysis involving the fitting of SARIMA models with precipitation, evaporation, and natural runoff as regressors to explain variations in Lake Malawi’s water level. The approach began with modeling without differencing. After observing suboptimal results, differencing was introduced to stabilize the series and improve model reliability. Model selection prioritized the lowest AIC values, with inverse roots plots and residual diagnostics used to assess stability and stationarity.

Initially, we fit SARIMA models without differencing. A grid search was conducted across various \((p, q)\) pairs and seasonal orders of (1,0,1), (0,0,1), and (1,0,0). Among these, the seasonal order of (1,0,1) yielded the lowest AIC values. Consequently, we focus on this seasonal order and present its corresponding AIC table below.

stats::arima()-Produced AIC Table for Signal Plus Noise Model
	MA0	MA1	MA2	MA3	MA4
AR0	4,310.70	4,019.08	3,835.52	3,715.60	3,637.01
AR1	3,523.31	3,491.39	3,484.51	3,486.33	3,488.32
AR2	3,484.10	3,485.85	3,486.38	3,488.45	3,490.34
AR3	3,485.76	3,486.74	3,488.42	3,490.28	3,492.34
AR4	3,486.09	3,487.70	3,489.85	3,491.79	3,488.61
^a All models contain a $(1, 0, 1)_{12}$ seasonal term.

In this table, we observe that the best-performing model is SARMA(2,0) × (1,1)₁₂ with the lowest AIC value of 3484.076, followed closely by SARMA(1,2) × (1,1)₁₂ with an AIC of 3484.506. Generally, adding a single parameter should increase the AIC by no more than two units and larger deviations may indicate potential maximization or evaluation errors during model fitting. However, certain transitions in the table deviate from this expectation. Despite these irregularities, the AIC values around the (2,0) region remain consistently low, reinforcing the choice of SARMA(2,0) × (1,1)₁₂ as the best model.

Despite the favorable AIC, we observe that the inverse roots plot for the SARMA(2,0) × (1,1)₁₂ model showed roots lying on the unit circle, indicating potential for non-stationarity. As discussed above, inverse roots plot is a diagnostic tool used to assess the stationarity and invertibility of ARMA models. For a model to be stationary and invertible, all inverse roots should lie inside the unit circle. Roots on or near the boundary suggest the model may not adequately capture the underlying process, potentially leading to unreliable forecasts.

Differencing

In our analysis, recognizing the non-stationarity present in the initial models, we applied differencing to stabilize the mean and variance, thereby enhancing the interpretability of the regression coefficients. This approach aligns with the findings of Kamalov (2021), who advocate for differencing as a means to achieve stationarity in time series data, ensuring that the models accurately reflect the underlying hydrological processes. In modeling Lake Malawi’s water level, differencing is justified for the following reasons:

• Long-term Trend Removal: The lake’s water levels show persistent trends from climatic or land-use changes. First differencing removes these long-term shifts, allowing regression coefficients to capture relevant short-term variations rather than spurious relationships.

• Seasonal Effects: Distinct wet and dry seasons influence water levels annually. Seasonal differencing (lag of 12 months) accounts for these predictable cycles, enabling the model to focus on irregular fluctuations and the direct effects of precipitation, evaporation, and runoff.

• In hydrological studies, according to Tunnicliffe Wilson (2016), failing to difference when data exhibit trends and seasonality can lead to models that overstate or misrepresent the influence of variables like precipitation or evaporation. For Lake Malawi, which is sensitive to seasonal rains and evaporation cycles, proper differencing ensures that regression results reflect real-time relationships rather than long-term shifts or seasonal patterns.

The plot of the first and seasonal differenced water level data shows fluctuations centered around a stable mean with no visible trend, indicating that differencing successfully removed non-stationarity. The variability appears consistent over time, suggesting the variance has been stabilized. Unlike the original series, this differenced series focuses on short-term changes. The absence of persistent upward or downward movements and the uniform spread of values justify the use of both first and seasonal differencing, making the data suitable for SARIMA modeling.

stats::arima()-Produced AIC Table for Signal Plus Noise Model
	MA0	MA1	MA2	MA3	MA4
AR0	3,342.51	3,321.35	3,317.36	3,319.36	3,321.35
AR1	3,316.78	3,318.63	NA	3,321.36	3,323.36
AR2	3,318.54	3,319.66	3,321.35	3,323.19	3,325.27
AR3	3,319.23	3,321.23	3,323.06	3,325.05	3,327.04
AR4	3,321.23	3,323.23	3,321.75	3,324.02	3,327.85
^a All models contain a $(1, 1, 1)_{12}$ seasonal term.
^b NA values correspond to models with failed convergence due to non-finite difference value.

Using differencing, a second grid search was conducted using seasonal order (1,1,1). The SARIMA(1,1,0) × (1,1,1)₁₂ model had the lowest AIC value of 3,316.78 and was selected as the final model.

arima_110x111 <- arima(data$W_level_mm, 
                       order = c(1,1,0), 
                       seasonal = list(order = c(1,1,1), period = 12), 
                       xreg = data[, c("Precip", "Evap", "Nat_runof")], 
                       method = "ML", optim.control = list(maxit = 1000))

The inverse roots plot for the fitted SARIMA(1,1,0) × (1,1,1)₁₂ model shows a single inverse AR root well inside the unit circle, which indicates that the model is stationary after applying both first and seasonal differencing. This confirms that the differencing effectively addressed non-stationarity and that the AR component is stable. The absence of roots near or on the boundary suggests that the model should produce reliable forecasts without explosive behavior or persistent trends. This observation aligns with the findings of Ikughur, Tersoo, and Oyewole (2015), which emphasized the importance of residual analysis in time series modeling and forecasting. Their study demonstrated that models with residuals exhibiting white noise characteristics—indicative of uncorrelated and homoscedastic residuals—tend to provide more accurate and reliable forecasts.

The diagnostic plots were used to assess the adequacy of the fitted SARIMA(1,1,0) × (1,1,1)₁₂ model.

The residual time series plot (left) shows residuals fluctuating around zero, indicating no obvious trend or seasonality. However, there are some large spikes early in the series, suggesting possible outliers or periods of volatility.

The Normal Q-Q plot (right) checks the normality of residuals. While most points lie along the 45-degree line, deviations occur at the tails, indicating mild departures from normality. Such deviations are common but should be noted when making probabilistic forecasts.

The ACF plot of residuals shows that most autocorrelations lie within the significance bounds. The absence of significant spikes suggests the model has adequately captured the data’s temporal dependence structure.

Considering the combination of stable residuals, minimal autocorrelation, and reasonably normal error distribution, the SARIMA(1,1,0) × (1,1,1)₁₂ model provides a good balance of interpretability and predictive performance. Despite minor imperfections, these diagnostics indicate no major issues that would necessitate a more complex model, justifying its selection as the final model.

Coefficients

Model Coefficients
Feature	Coefficient	Standard Error
Precipitation	−0.0161	0.0980
Evaporation	−0.9463	0.3697
Natural Runoff	−0.2377	0.1349

In the table above, we observe all coefficients are negative - i.e., increases in precipitation, evaporation, and natural runoff for a given month are associated with decreased water levels that month. These results were initially surprising. While we hypothesized evaporation and runoff would be associated with decreased water levels (given they involve water leaving a system), precipitation would seem to logically be connected with increased water levels.

Upon further examination, we see a relatively large standard error for the effect of precipitation (compared with its coefficient). This raised the question of whether precipitation was a material predictor of water level (regardless of sign). Wanting a hypothesis test that reflected profile likelihood-based CIs, we used a likelihood ratio test for the following: \[H_0: \beta_{precipitation} = 0,\text{ } H_0: \beta_{precipitation} \neq 0.\]

The difference in log-likelihoods between the two models was \(\approx 0.014\), inducing a p-value of \(\approx 0.869\). Subsequently, we would not reject the null hypothesis that precipitation does not have a linear association with water level.

As a sanity check, we also examined whether or not all three variables were associated with water level (jointly testing \(\beta_{precipitation}, \beta_{evaporation}, \beta_{runoff}=0\)). We did so via likelihood ratio test, as above, and obtained a p-value of \(\approx 0.003\). Therefore, we would conclude that these three variables together provide explanatory value.¹¹

Conclusions

After attempting to forecast evaporation, precipitation, and natural runoff in the Lake Malawi region, we found that precipitation was the most difficult variable of the three to forecast. However, quite surprisingly, we found that precipitation had little association with water level. Weighting the MAE for each variable by its impact on water level, reveals the following.

Mean Absolute Error Weighted by Impact (mm)

Variable	Raw MAE	Weighted MAE
Evaporation	6.32	5.98
Precipitation	32.27	0.52
Runoff	15.66	3.72

While evaporation was the easiest to predict and precipitation the most difficult, weighting by their influence on water level shows that forecasting errors for evaporation pose a more significant challenge to predicting water levels than those for precipitation.

Again, we found this quite surprising. However, further research provided context and intuition for the observed patterns. As noted in (Mtilatila, Bronstert, and Vormoor 2022), the Malawian government installed a barrage in the 1960s to regulate outflow for agricultural and hydroelectric purposes, directly influencing lake levels. This artificial intervention likely disrupted the natural relationship between precipitation and water level; for example, following heavy rainfall, the government may release more water from the lake, weakening the direct link between precipitation and water level.¹² Given this regulatory context, we caution against extrapolating our findings beyond Lake Malawi.

More broadly, our findings underscore the importance of considering human interventions when modeling hydrological systems. While traditional forecasts often emphasize precipitation as a key driver of water levels, our results highlight the need for a more nuanced approach—one that accounts for both environmental processes and regulatory controls. Future research should further explore how such interventions shape hydrological dynamics, particularly in regions where water management plays a critical role in resource allocation.

Scholarship: Comparison with Past Projects and Lessons from Peer Review

Our project builds on previous DATASCI/STATS 531/631 projects, specifically

1. Precipitation in Detroit (Winter 2024)¹³

2. CO2 Emissions in the US (Winter 2024)¹⁴

Like (1) above, we explore time series models for hydrological systems. However, this project focused on establishing a formal, data-driven, approach to constructing a seasonal model (e.g., exploring a smoothed periodogram to investigate plausible periods of seasonality). In contrast, as discussed above, we chose to rely on scientific literature/domain knowledge and assume yearly fluctuations (furthermore, this pattern was clear in the data). Furthermore, we explored several different aspects of the local water cycle, whereas (1) primarily discussed precipitation. Additionally, we found log transforming our precipitation data to be useful in model building, whereas this was not the case in (1).

Like (2) above, we consider data partitions for our analysis. However, the motivation behind the data splitting in (2) appears to be analysis of persistence of a trend. I.e., they claim to observe different trends in recent years compared with the initial years in their data; subsequently, they investigate the sensitivity of their findings to the choice of “start year” for their training data. On the other hand, as discussed in the Model Evaluation section, we partition the data to achieve a test/train split. This choice allows us to use our observed MAE values as valid estimators of true model accuracy. Furthermore, as suggested in the peer review for (2), we chose to investigate different seasonal \((P,Q)_{12}\) structures to help address residual heteroscedasticity.

Beyond these two projects in particular, we find similarities with our project and most past projects in our focus on ARMA/SARMA modeling, with AIC-based approaches for model selection. Our project mainly deviates from most others in the following senses:

• Motivation for the research question (looking at components of the variable of interest to investigate potential forecasting obstacles).

• Specific focus on how optim() maximizes the likelihood in model fitting, and consequences of modifying default approaches.

• Explicit analysis of forecasting metrics of our models.

References

ChatGPT was used to help polish individual sentences/paragraphs.

Bhave, Ajay G., Lauren Bulcock, Suraje Dessai, Declan Conway, Graham Jewitt, Andrew J. Dougill, Seshagiri Rao Kolusu, and David Mkwambisi. 2020. “Lake Malawi’s Threshold Behaviour: A Stakeholder-Informed Model to Simulate Sensitivity to Climate Change.” Journal of Hydrology 584 (May): 124671. https://doi.org/10.1016/j.jhydrol.2020.124671.

Bootsma, Harvey A., and Robert E. Hecky. 2003. “A Comparative Introduction to the Biology and Limnology of the African Great Lakes.” Journal of Great Lakes Research 29 (January): 3–18. https://doi.org/10.1016/s0380-1330(03)70535-8.

Bootsma, Harvey A., and Sven Erik Jørgensen. 2013. “Lake Malawi/Nyasa: Experience and Lessons Learned Brief.” In. https://api.semanticscholar.org/CorpusID:133288245.

Chiwaula, Levison, Daniel Jamu, R. Chaweza, and Joseph Nagoli. 2012. “The Structure and Margins of the Lake Chilwa Fisheries in Malawi: A Value Chain Analysis,” January.

Chou, Chia, John C. H. Chiang, Chia-Wei Lan, Chia-Hui Chung, Yi-Chun Liao, and Chia-Jung Lee. 2013. “Increase in the Range Between Wet and Dry Season Precipitation.” Nature Geoscience 6 (4): 263–67. https://doi.org/10.1038/ngeo1744.

Hyndman, Rob. “Transforming Data with Zeros.” https://robjhyndman.com/hyndsight/transformations/.

Ikughur, Atsua, Uba Tersoo, and Adeniyi Oyewole. 2015. “Application of Residual Analysis in Time Series Model Selection Application of Residual Analysis in Time Series Model Selection.” Journal of Economics and Econometrics 4 (January): 41–53.

Jury, M. R., and N. D. Mwafulirwa. 2002. “Climate Variability in Malawi, Part 1: Dry Summers, Statistical Associations and Predictability.” International Journal of Climatology 22 (11): 1289–1302. https://doi.org/10.1002/joc.771.

Jury, Mark R. 2014. “Malawi’s Shire River Fluctuations and Climate.” Journal of Hydrometeorology 15 (5): 2039–49. https://doi.org/10.1175/jhm-d-13-0195.1.

Kamalov, Firuz. 2021. “A Note on Time Series Differencing.” Gulf Journal of Mathematics 10 (May): 50–56. https://doi.org/10.56947/gjom.v10i2.609.

Kumambala, Patsani G. 2010. “Water Balance Model of Lake Malawi and Its Sensitivity to Climate Change.” The Open Hydrology Journal 4 (1): 152–62. https://doi.org/10.2174/1874378101004010152.

Lancaster, N. 1979. “The Physical Environment of Lake Chilwa.” In Lake Chilwa: Studies of Change in a Tropical Ecosystem, 17–40. Springer.

Lyons, Robert P., Christopher A. Scholz, Andrew S. Cohen, John W. King, Erik T. Brown, Sarah J. Ivory, Thomas C. Johnson, et al. 2015. “Continuous 1.3-Million-Year Record of East African Hydroclimate, and Implications for Patterns of Evolution and Biodiversity.” Proceedings of the National Academy of Sciences 112 (51): 15568–73. https://doi.org/10.1073/pnas.1512864112.

Mbaye, Mamadou Lamine, Stefan Hagemann, Andreas Haensler, Tobias Stacke, Amadou Thierno Gaye, and Abel Afouda. 2015. “Assessment of Climate Change Impact on Water Resources in the Upper Senegal Basin (West Africa).” American Journal of Climate Change 04 (01): 77–93. https://doi.org/10.4236/ajcc.2015.41008.

Mtilatila, Lucy, Axel Bronstert, Pallav Shrestha, Peter Kadewere, and Klaus Vormoor. 2020. “Susceptibility of Water Resources and Hydropower Production to Climate Change in the Tropics: The Case of Lake Malawi and Shire River Basins, SE Africa.” Hydrology 7 (3): 54. https://doi.org/10.3390/hydrology7030054.

Mtilatila, Lucy, Axel Bronstert, and Klaus Vormoor. 2022. “Temporal Evaluation and Projections of Meteorological Droughts in the Greater Lake Malawi Basin, Southeast Africa.” Zweitveröffentlichungen Der Universität Potsdam : Mathematisch-Naturwissenschaftliche Reihe; 1287. https://doi.org/10.25932/PUBLISHUP-57128.

Nicholson, Sharon E. 2017. “Climate and Climatic Variability of Rainfall over Eastern Africa.” Reviews of Geophysics 55 (3): 590–635. https://doi.org/10.1002/2016rg000544.

Njaya, Friday, Katherine A. Snyder, Daniel Jamu, John Wilson, Clive Howard-Williams, Edward H. Allison, and Neil L. Andrew. 2011. “The Natural History and Fisheries Ecology of Lake Chilwa, Southern Malawi.” Journal of Great Lakes Research 37 (January): 15–25. https://doi.org/10.1016/j.jglr.2010.09.008.

Tunnicliffe Wilson, Granville. 2016. “Time Series Analysis: Forecasting and Control,5th Edition, by George e. P. Box, Gwilym m. Jenkins, Gregory c. Reinsel and Greta m. Ljung, 2015. Published by John Wiley and Sons Inc., Hoboken, New Jersey, Pp. 712. ISBN: 978-1-118-67502-1.” Journal of Time Series Analysis 37 (March): n/a–. https://doi.org/10.1111/jtsa.12194.

Wheeler, Jesse, Noel McAllister, and Dhajanae Sylvertooth. 2024. Arima2: Likelihood Based Inference for ARIMA Modeling. https://CRAN.R-project.org/package=arima2.

---

1. https://www.climateengine.org/↩︎

2. Evaporation, precipitation, and runoff were obtained via Climate Engine. However, water level (not publicly available) was obtained via the Malawi government. We thank
   Mexford Chimwemwe Mulumpwa
   Department of Aquaculture and Fisheries
   Lilongwe University of Agriculture and Natural Resources
   P.O. Box 219
   Lilongwe 3
   for the water level data.↩︎

3. Note: We could have examined a (smoothed) periodogram to confirm that the dominant oscillations took the form of annual cycles. However, the data shows such strong annual seasonality (backed by scientific/cultural consensus) that we did not feel a data-driven exploration would be valuable here.↩︎

4. On the advice of Jesse Wheeler.↩︎

5. See https://www.rdocumentation.org/packages/arima2/versions/3.3.0/topics/arima↩︎

6. For the purposes of AIC tables, we used the above method as the initial approach. If this failed, then we also set method="ML", which uses a maximum likelihood fitting method (rather than using conditional-sum-of-squares to find starting values for the likelihood maximization).↩︎

7. Full table omitted for brevity.↩︎

8. This model did not have any AR/MA roots on/in the unit circle.↩︎

9. The models \(SARMA(2,4)\times(1,1)_{12}\), \(SARMA(3,1)\times(1,1)_{12}\), and \(SARMA(2,1)\times(1,1)_{12}\) had lower AIC values; however, all had unit roots in the AR and/or MA polynomials.↩︎

10. I.e., we invert the transformations.↩︎

11. Note that we would still reject the null here when using a Bonferroni-adjusted testing threshold (\(\alpha'=0.025\)) to account for having run two hypothesis tests above.↩︎

12. It may be difficult to quantify this exactly, as the government has not always followed their stated outflow policy (Mtilatila, Bronstert, and Vormoor 2022).↩︎

13. https://ionides.github.io/531w24/midterm_project/project05/blinded.html↩︎

14. https://ionides.github.io/531w24/midterm_project/project06/blinded.html↩︎