Javascript is required
Anderson, C. & Ryan, L. (2017). A comparison of spatio-temporal disease mapping approaches including an application to ischaemic heart disease in New South Wales, Australia. Int. J. Environ. Res. Public Health, 14(2), 146. [Google Scholar] [Crossref]
FAO. (2023). The State of Food Security and Nutrition in the World 2023. Rome, Italy: FAO; IFAD; UNICEF; WFP; WHO. [Google Scholar] [Crossref]
Hidayat, Y., Purwandari, T., Ratnasari, D., Sukono, Saputra, J., & Subiyanto. (2021). Objective method for determining the importance of unprecedented restlessness as a rice crisis indicator at the national level. Agronomy, 11(6), 1195. [Google Scholar] [Crossref]
Hyndman, R. J. & Athanasopoulos, G. (2018). Forecasting: Principles and Practice (2nd edition). OTexts: Melbourne, Australia. https://otexts.com/fpp2/ [Google Scholar]
International Research Institute for Climate & Society. (2016). Map of El Nino. https://iridl.ldeo.columbia.edu/maproom/IFRC/FIC/elninorain.html#tabs-2 [Google Scholar]
Kano, K., Kano, T., & Takechi, K. (2022). The price of distance: Pricing-to-market and geographic barriers. J. Econ. Geogr., 22(4), 873–899. [Google Scholar] [Crossref]
Kim, S. (2024). Transmission of economic growth between neighboring countries in Europe. Appl. Econ., 1–16. [Google Scholar] [Crossref]
LeSage, J. P. (2008). An introduction to spatial econometrics. Rev. Econ. Ind., 19–44. [Google Scholar] [Crossref]
Ljung, G. M. & Box, G. E. P. (1978). On a measure of lack of fit in time series models. Biometrika, 65(2), 297–303. [Google Scholar] [Crossref]
Lu, G. Y. & Wong, D. W. (2008). An adaptive inverse-distance weighting spatial interpolation technique. Comput. Geosci., 34(9), 1044–1055. [Google Scholar] [Crossref]
Makridakis, S., Wheelwright, S. C., & Hyndman, R. J. (1998). Forecasting: Methods and Applications (Third Edition). New York: John Wiley & Sons. [Google Scholar]
Mardia, K. V. (1970). Measures of multivariate skewness and kurtosis with applications. Biometrika, 57(3), 519–530. [Google Scholar] [Crossref]
Mustafa, M., Yuyun, H., & Ismail, M. (2015). Developing unprecedented restlessness as the new strong indicator of rice crisis at national level. Glob. J. Pure Appl. Math., 11(6), 4139–4160. [Google Scholar]
Oliver, M. A. & Webster, R. (2015). Basic Steps in Geostatistics: The Variogram and Kriging. Springer. [Google Scholar] [Crossref]
Pfeifer, P. E. & Deutsch, S. J. (1980). A three-stage iterative procedure for space-time modeling. Technometrics, 22(1), 35–47. [Google Scholar] [Crossref]
Rahmawati, S. (2006). Status perkembangan perbaikan sifat genetik padi menggunakan transformasi argobacterium. J. Agrobiogen, 2(1), 36–44. [Google Scholar]
Ruchjana, B. N. (2019). Pengembangan model spatio temporal dan aplikasinya. In Prosiding Seminar Nasional Matematika Dan Pendidikan Matematika, UIN Raden Intan Lampung, 2(1). https://ejournal.radenintan.ac.id/index.php/pspm/article/view/4295 [Google Scholar]
Scaife, A., Guilyardi, E., Cain, M., & Gilbert, A. (2019). What is the El Niño–Southern Oscillation? Weather, 74(7), 250–251. [Google Scholar] [Crossref]
Soni, A. & Singh, P. (2022). A review paper on El Nino. Int. J. Innov. Res. Eng. Manag., 9(1), 269–272. [Google Scholar] [Crossref]
Tobler, W. R. (1970). A computer movie simulating urban growth in the Detroit region. Econ. Geogr., 46, 234–240. [Google Scholar] [Crossref]
Wong, D. W. S. & Lee, J. (2005). Statistical Analysis of Geographic Information with ArcView GIS and ArcGIS. John Wiley & Sons, Inc. [Google Scholar]
Search
Open Access
Research article

Forecasting Rice Crisis Dynamics in Southeast Asia Using a Spatio-Temporal Autoregressive Model Based on the Restlessness Indicator

yuyun hidayat,
bertho tantular,
cosmas david ananda manurung*
Department of Statistics, Padjadjaran University, 45363 Jatinangor, Indonesia
Organic Farming
|
Volume 11, Issue 1, 2025
|
Pages 13-38
Received: 02-06-2025,
Revised: 03-13-2025,
Accepted: 03-20-2025,
Available online: 03-30-2025
View Full Article|Download PDF

Abstract:

The rice crisis represents a significant threat to food security and economic stability in Southeast Asia, a region where rice serves as the primary staple for the majority of the population. This crisis is exacerbated by a confluence of factors, including climate change, crop failures, and restrictive export policies, as exemplified by the El Niño phenomenon and India’s 2023 rice export ban. Rising rice prices have been linked to increased social unrest, with the potential to trigger widespread demonstrations across affected nations. To proactively address this issue, the restlessness indicator was introduced as a predictive tool, integrating key variables such as rice prices, consumption patterns, and per capita income. This study employs a Spatio-Temporal Autoregressive (STAR) model to forecast restlessness values across six Southeast Asian countries—Indonesia, the Philippines, Thailand, Vietnam, Malaysia, and Cambodia—from 2024 to 2028. The STAR (5,1) model was identified as the optimal framework, achieving a Mean Absolute Percentage Error (MAPE) of 15.1%. The forecasting results indicate that none of the analyzed countries are projected to enter a state of unprecedented restlessness during the specified period, suggesting that no severe rice crisis is anticipated within this timeframe. These findings provide critical insights for policymakers and stakeholders, enabling the development of preemptive strategies to mitigate potential food security challenges. The study underscores the utility of the restlessness indicator as a robust tool for monitoring and forecasting rice-related crises, contributing to the broader discourse on sustainable food systems in Southeast Asia.

Keywords: Rice crisis, Spatio-temporal autoregressive model, Food security, Southeast Asia, Restlessness indicator, Forecasting

1. Introduction

1.1 Background

The majority of the world's population, especially in developing countries such as Indonesia, relies on rice as their primary food source to meet daily nutritional needs (R​a​h​m​a​w​a​t​i​,​ ​2​0​0​6). According to the data from the Foreign Agricultural Service (FAS) of the United States Department of Agriculture (USDA), rice consumption in Asia accounts for 80% of the total global rice consumption. Figure 1 shows the comparison of rice consumption between Asia and the world.

Figure 1. Comparison of rice consumption between Asia and the world

In 2023, most regions in Asia experienced dry conditions or droughts, accompanied by a drastic reduction in rainfall. This phenomenon, known as El Niño, is a major climate disturbance that occurs every two to eight years in the equatorial Pacific Ocean (S​o​n​i​ ​&​ ​S​i​n​g​h​,​ ​2​0​2​2). The impact of El Niño results in increased rainfall in Latin America (S​c​a​i​f​e​ ​e​t​ ​a​l​.​,​ ​2​0​1​9), while the opposite occurs in Australia and Asia. Although not all parts of Asia experience drought, some regions receive excessive rainfall due to El Niño, such as northern and southwestern India.

Extreme climate change contributes to crop failures, notably in India. High rainfall in India has caused flooding in several regions, including major rice-producing areas. According to reports from India's Ministry of Agriculture, 49.8 million hectares of land (approximately 15.2% of India's total area) were flooded, with a significant portion being rice fields. Figure 2 shows the rainfall pattern related to El Niño across various areas worldwide (I​n​t​e​r​n​a​t​i​o​n​a​l​ ​R​e​s​e​a​r​c​h​ ​I​n​s​t​i​t​u​t​e​ ​f​o​r​ ​C​l​i​m​a​t​e​ ​&​ ​S​o​c​i​e​t​y​,​ ​2​0​1​6).

Figure 2. Rainfall pattern related to El Niño across various areas worldwide
Source: I​n​t​e​r​n​a​t​i​o​n​a​l​ ​R​e​s​e​a​r​c​h​ ​I​n​s​t​i​t​u​t​e​ ​f​o​r​ ​C​l​i​m​a​t​e​ ​&​ ​S​o​c​i​e​t​y​ ​(​2​0​1​6​)

Due to crop failures, India officially imposed a rice export ban starting from July 20, 2023, leading to a surge in rice prices (inflation) across Asia. According to Consumer News and Business Channel (CNBC), the price of rice in Indonesia as of September 2023 reached IDR 14,290/kg, marking a 14.5% increase from IDR 12,480/kg in September of the previous year. Based on observations from Kompas, the number of regions experiencing rice price increases rose from 300 districts/cities in the first week of September 2023 to 341 districts/cities in the second week of September 2023. According to the Voice of America (VOA), the President of the Philippines implemented price controls on September 5, 2023, after rice prices reached their highest point in 14 years. According to Nation Thailand (2023), in Thailand, rice prices reached a 17-year high in September 2023, prompting the government to urge farmers to plant only one rice crop to conserve water amid limited rainfall.

A rice crisis can have far-reaching consequences across various levels of society, particularly affecting impoverished populations, government policies, and global trade dynamics. For impoverished populations, rice is often a staple food and a primary source of calories. A shortage or sharp increase in rice prices can lead to food insecurity, malnutrition, and heightened poverty levels. Families may be forced to allocate a larger portion of their income to purchasing rice, leaving less for other essential needs such as healthcare and education. This situation can exacerbate inequality and deepen the vulnerability of already marginalized groups (F​A​O​,​ ​2​0​2​3). According to BBC News in February 2024, many low-income households in Indonesia opted to buy subsidized rice (SPHP) due to high rice prices. People queued from early morning to obtain subsidized rice distributed by Bulog, which ran out within an hour, leaving many without access due to limited supply. In some areas, frustration and dissatisfaction over rising rice prices led to protests. According to CNBC Indonesia, workers staged a demonstration on Jalan Merdeka Barat on September 25, 2023, due to continuously rising rice prices. Students in Sukabumi also protested in front of the DPRD building in March 2024 over skyrocketing rice prices in the region.

At the governmental level, a rice crisis often prompts immediate policy responses to stabilize supply and prices. Governments may implement export bans, subsidies, or price controls to protect domestic consumers. However, such measures can sometimes backfire, leading to market distortions, reduced incentives for farmers, and strained international trade relations. For instance, export restrictions by major rice-producing countries can disrupt global supply chains, further inflating prices and creating tensions among trading partners. Additionally, governments may need to increase spending on social safety nets, such as food aid or cash transfers, which can strain national budgets and divert resources from other critical areas like infrastructure or education.

On a global scale, a rice crisis can disrupt international trade and geopolitical stability. The crisis can intensify competition for resources, leading to potential conflicts or alliances formed around securing food supplies. Thus, the ripple effects of a rice crisis extend beyond immediate food shortages, influencing economic, political, and environmental systems worldwide.

Climate change, crop failures, and export restrictions can make rice less accessible or cause shortages in rice-consuming countries like Indonesia. If left unchecked, this situation could lead to a rice crisis. The unpredictable factors that can harm nations and communities can be mitigated if governments anticipate potential scenarios and implement appropriate policies. Forecasting rice conditions in each country (especially in Asia) can serve as a consideration for the Ministry of Foreign Affairs of the Republic of Indonesia in formulating policies. These policies could include exporting rice to countries in need to generate economic benefits or imposing import or export restrictions to ensure domestic rice availability. Given the aforementioned factors, forecasting rice conditions using the restlessness indicator as a strong predictor of a national rice crisis is essential for policy formulation.

Restlessness (R) is an indicator that reflects a country's rice situation, derived from rice prices, rice consumption, and per capita income (H​i​d​a​y​a​t​ ​e​t​ ​a​l​.​,​ ​2​0​2​1). An increase in the restlessness value in a given year signifies greater difficulty in securing rice during that period.

Figure 3. Restlessness value in Southeast Asia in 2011

Based on Figure 3, the restlessness value in Southeast Asia in 2011 suggests a spatial effect, where neighboring countries tend to exhibit similar restlessness values. For example, the restlessness values of countries bordering Cambodia are similar (indicated by dark colors). Thailand and the Philippines, which are close to Malaysia, also show similar restlessness values to Malaysia. This phenomenon aligns with Tobler's First Law, which states that "everything is related to everything else, but near things are more related than distant things."

Figure 4 shows the temporal trend of restlessness values in Southeast Asian countries from 1992 to 2023. Indonesia, Brunei, the Philippines, Malaysia, and Thailand exhibit similar and stationary trends, while Cambodia shows fluctuations between 1995 and 2000. These fluctuations in Cambodia result from the country's economic and political dynamics. The per capita income used to calculate restlessness in Cambodia dropped significantly between 1995 and 2000, leading to a spike in restlessness values. Figure 5 shows the temporal trend of restlessness values in Southeast Asian countries from 1992 to 2023 (excluding Cambodia and Vietnam).

Given the suspected spatial and temporal effects in restlessness values, a forecasting method that accounts for both spatial and temporal dependencies is required. Therefore, this study forecasts restlessness values across multiple locations using the STAR model, a spatio-temporal approach. The forecast restlessness values were used to determine unprecedented restlessness in this study, which can indicate a potential rice crisis.

Figure 4. Restlessness values in Southeast Asia (1992-2023)
Figure 5. Restlessness values in Southeast Asian countries (1992-2023, excluding Cambodia and Vietnam)
1.2 Research Assumptions

The assumptions of this study state that the autoregressive coefficient of the STAR model remains constant for each location, meaning that the temporal pattern in Indonesia, the Philippines, Thailand, Vietnam, Malaysia, and Cambodia is considered homogeneous.

2. Methodology

2.1 Data Source

The data used in this study is secondary data sourced from the websites of the Food and Agriculture Organization (FAO) of the United Nations, the FAS, and the World Bank. This data consists of annual time-series data for six Southeast Asian countries (Indonesia, Philippines, Thailand, Vietnam, Malaysia, and Cambodia) spanning 32 years, from 1992 to 2023. The key variables include:

a) Producer rice price

The producer price for rice refers to the price received by farmers for their primary harvest collected at the initial selling point (the price paid at the farm gate). This data is expressed in $US/kg, with each year’s value representing the average price received by farmers throughout that year. Producer rice price data can be accessed through the link: (https://www.fao.org/faostat/en/#data/PP).

b) Income per capita

Income per capita, expressed in $US/person, represents the average income earned by an individual in a country over a specific period. This value is calculated by dividing the total national income (or GDP) by the total population. Income per capita data can be accessed through the link: (https://data.worldbank.org/indicator/NY.GDP.PCAP.CD).

c) Rice consumption per capita

Rice consumption per capita, expressed in kg/year, indicates the average amount of rice consumed per person throughout the year in a given country. Rice consumption per capita data can be accessed through the link: (https://apps.fas.usda.gov/psdonline/app/index.html#/app/home).

2.2 Restlessness

The variable used in this study is restlessness. Restlessness is an indicator used to predict future rice crises. It can be understood as unease, anxiety, or uncertainty in the context of rice availability. The significant increase in rice prices from year to year characterizes a significant increase in restlessness. Restlessness is obtained by dividing the rice price (USD/kg) by income per capita (USD/person) and then multiplying it by the rice consumption of each country. Meanwhile, unprecedented restlessness serves as a strong indicator of a rice crisis. The presence of unprecedented restlessness indicates an impending rice crisis in a country, and its absence suggests stability in rice availability. The formula for calculating restlessness is as follows (H​i​d​a​y​a​t​ ​e​t​ ​a​l​.​,​ ​2​0​2​1):

$R=\frac{P}{I} * C$
(1)

where, R is the restlessness value, P is the producer rice price ($US/kg), I is the income per capita ($US/person), and C is the rice consumption (kg/year).

This function is based on the idea that a rice crisis cannot be observed solely from price movements. High rice prices in a given year do not necessarily indicate a crisis if the per capita income also increases proportionally, as people can still afford their basic food needs. However, if rice prices rise while per capita income remains stagnant or decreases, the burden on consumers increases, leading to heightened restlessness.

Additionally, the consumption variable (C) accounts for population growth trends. If rice consumption is low in a particular year, the level of anxiety regarding rice availability is also reduced (M​u​s​t​a​f​a​ ​e​t​ ​a​l​.​,​ ​2​0​1​5). Therefore, restlessness serves as an early warning indicator for potential rice crises, providing policymakers with insights to take preventive measures before unrest occurs. Unprecedented restlessness is determined using a control chart, where values exceeding the Upper Control Limit (UCL) indicate a crisis scenario.

2.3 Data Preprocessing

Missing data points were addressed using spatial interpolation via the Inverse Distance Weighting (IDW). The IDW method estimates missing values based on data similarity between neighboring countries. The principle of IDW is that points closer to the missing data location have a greater influence than those farther away. The IDW formula (L​u​ ​&​ ​W​o​n​g​,​ ​2​0​0​8) is given as follows:

$y_0=\sum_{i=1}^n \lambda_i y_i$
(2)
$\lambda_i=\frac{\frac{1}{d_{0 i}^\alpha}}{\sum_i^n \frac{1}{d_{0 i}^\alpha}}$
(3)

where, $y_0$ is the estimated value at the missing data location, $y_i$ is the observed value at location $I, \lambda_i$ is the weight assigned to location $I, d_{0 i}^\alpha$ is the geographical distance between location 0 (missing data location) and location $i$ with the power of $\alpha$.

The $\alpha$ parameter is specified as a geometric form for the weight while other specifications are possible. This specification implies that if $\alpha$ is larger than 1, the so-called distance-decay effect will be more than proportional to an increase in distance, and vice versa. Thus, small $\alpha$ tends to yield estimated values as averages of $y_i$ in the neighborhood, while large $\alpha$ tends to give larger weights to the nearest points and increasingly down-weights points farther away. In this study, $\alpha$ = 1 was used to balance the weighting effect between nearby and farther locations.

2.4 Data Analysis
2.4.1 STAR

A spatio-temporal model is a model that integrates spatial and temporal dependencies into a time-series dataset. The STAR model is a statistical model used for forecasting spatio-temporal time series data, developed by Pfeifer and Deutsch in 1979. This model is designed to capture the relationship between time-series data recorded simultaneously across multiple locations while accounting for both spatial and temporal dependencies (R​u​c​h​j​a​n​a​,​ ​2​0​1​9). The STAR model is an extension of the Vector Autoregressive (VAR) model, which only considers temporal dependencies. The STAR (p, λk) model is expressed by the following equation (P​f​e​i​f​e​r​ ​&​ ​D​e​u​t​s​c​h​,​ ​1​9​8​0):

$z(t)=\sum_{k=1}^p \sum_{l=0}^{\lambda_k}\left[\phi_{k l} W^{(l)} z(t-k)\right]+\epsilon(t)$
(4)

where, $z(t)$ is the vector $(N \times 1)$ representing $N$ locations at time $t, z(t-p)$ is the vector $(N \times 1)$ representing $N$ locations at the time $(t-p), \lambda_k$ is the spatial order of the autoregressive model of order $k, \phi_{k l}$ is the STAR model parameter for time lag $k$ and spatial lag $l, W^{(l)}$ is the weight matrix of size $(N \times N)$ at spatial lag $l$, where the diagonal elements are zero, and each row sums to one; and $\epsilon_t$ is the error vector of size $(N \times 1)$ at time $t$.

The STAR model was chosen over alternative spatio-temporal models such as the Spatio-Temporal Moving Average (STMA) and Spatio-Temporal Autoregressive Moving Average (STARMA) models due to its superior ability to capture complex spatial and temporal dependencies. Unlike STMA, which primarily captures short-term fluctuations through moving averages, STAR explicitly incorporates autoregressive terms. This means that STAR accounts for how past observations influence future values, making it more suitable for forecasting long-term trends in restlessness values. This autoregressive nature is crucial in economic and agricultural contexts, where effects are often persistent and cumulative rather than transient (L​e​S​a​g​e​,​ ​2​0​0​8). Another advantage of STAR over STMA is its improved forecasting stability. The inclusion of autoregressive components ensures greater consistency in predictions compared to STMA, which can be more sensitive to sudden shocks and localized variations. By considering past values in both time and space, STAR provides smoother and more reliable forecasts, which is critical for policy formulation and crisis mitigation.

While STARMA models incorporate both autoregressive and moving average components in time and space, they tend to be computationally intensive due to the large number of parameters required for estimation (P​f​e​i​f​e​r​ ​&​ ​D​e​u​t​s​c​h​,​ ​1​9​8​0). STAR, on the other hand, achieves comparable predictive accuracy with fewer parameters, making it more practical for large-scale applications such as monitoring rice crises across multiple countries. The computational efficiency of STAR allows for quicker and more reliable model estimations, which is essential for timely decision-making.

The selection of STAR in this study is further supported by empirical evidence from the data, which exhibits significant partial autocorrelation in both temporal and spatial dimensions. The presence of nonzero lag effects in the Partial Autocorrelation Function (PACF) suggests that an autoregressive approach is more appropriate than a moving average model like STMA (H​y​n​d​m​a​n​ ​&​ ​A​t​h​a​n​a​s​o​p​o​u​l​o​s​,​ ​2​0​1​8). This characteristic justifies the use of STAR, as it effectively captures the persistence of restlessness values over time and across regions.

2.4.2 Spatial weights

Spatial weights are a crucial component used to capture spatial dependencies between locations in spatio-temporal modeling. Spatial dependencies explain how a variable at one location is influenced by the same variable at other locations. Two common types of spatial weights are contiguity-based and distance-based weights.

a) Contiguity-based weights: These weights are based on physical proximity or shared borders between spatial units, such as countries, provinces, or districts. This approach assumes that interactions are strongest among neighboring spatial units that share borders.

b) Distance-based weights: These weights use geographical distances between coordinates to determine the strength of spatial relationships. This method is more suitable when spatial units do not share borders but still exhibit interactions based on proximity.

In addition to weighting based on contiguity and distance, weighting can also use trade data-based weights, specifically inter-country rice trade data, as rice significantly influences restlessness itself. Trade data-based weights can be calculated using the following equation:

$w_{i j}=\frac{T_{i j}}{\sum_{j=1}^n T_{i j}}, i \neq j$
(5)

where, $w_{i j}$ is the spatial weight between locations $i$ and $j$, and $T_{i j}$ is the trade volume from location $i$ to location $j$.

However, there are data limitations for international trade data. Therefore, the use of weights with rice trade data is not possible. Therefore, in this study the most appropriate spatial weighting used is distance-based weighting.

The inverse distance weight is the most commonly used spatial weight based on distance. This weighting is applied based on the inverse of the geographical distance between locations. There is a single coordinate point representing the geographical center (centroid) of each research object. The weights can be calculated using the following equation:

$w_{i j}=\frac{\frac{1}{d_{i j}}}{\sum_{j=1}^n \frac{1}{d_{i j}}}, i \neq j$
(6)

where, $w_{i j}$ is the inverse distance weight between locations $i$ and $j$, and $d_{i j}$ is the geographical distance between locations $i$ and $j$.

The use of IDW is also supported by the strong correlation between distance and the import volume of several countries, such as the Philippines (correlation -0.6484), which imports more from closer countries, and Cambodia (correlation -0.8286). Additionally, distance is an important factor because rice movement is also influenced by transportation and logistics. Greater distances can lead to increased transportation time and costs, which in turn can affect rice prices for consumers (K​a​n​o​ ​e​t​ ​a​l​.​,​ ​2​0​2​2). Moreover, the restlessness value is also constructed from the income variable, which is closely related to the economic conditions of the country. Strong economic growth in a country can drive development and prosperity in neighboring countries through trade and investment relationships (K​i​m​,​ ​2​0​2​4).

2.4.3 Spatial autocorrelation testing (MoranST)

Spatial autocorrelation measures the extent to which a variable’s value at one location is related to the values at neighboring locations. Based on Tobler's First Law of Geography, "everything is related to everything else, but near things are more related than distant things" (T​o​b​l​e​r​,​ ​1​9​7​0). Spatial autocorrelation can be detected using the Moran’s I statistic, calculated as (A​n​d​e​r​s​o​n​ ​&​ ​R​y​a​n​,​ ​2​0​1​7):

${MoranST}=\frac{N T \sum_{i=1}^N \Sigma_{t=1}^T \Sigma_{j=1}^N \Sigma_{S=1}^T \tilde{w}_{(t t, j s)}\left(y_{i t}-\tilde{y}\right)\left(y_{j s}-\tilde{y}\right)}{\sum_{i=1}^N \Sigma_{t=1}^T \Sigma_{j=1}^N \Sigma_{S=1}^T \tilde{w}_{(i t, j s)} \sum_{i=1}^N \Sigma_{t=1}^T\left(y_{i t}-\tilde{y}\right)^2}$
(7)

where, $\tilde{y}$ is the average of observed variables from $N$ spatial units and $T$ temporal units, $\widetilde{w}_{(i t, j s)}$ is the elements of the weighting matrix for the spatio-temporal autocorrelation between $y_{i t}$ and $y_{j s}, y_{i t}$ is the value of observed variables of the $i$-th spatial unit and the $t$-th temporal unit, and $y_{j s}$ is the value of observed variables of the $j$-th spatial unit and the $s$-th temporal unit.

The hypothesis for Moran’s I test is:

H0: There is no spatial autocorrelation (random spatial distribution).

H1: There is spatial autocorrelation (non-random spatial distribution).

Moran’s I ranges from -1 to 1. If the Moran index value is in the range of 0 < I ≤ 1, this indicates positive spatial autocorrelation, meaning that adjacent locations tend to cluster spatially. Conversely, if the Moran index value is in the range of -1 ≤ I < 0, this indicates negative spatial autocorrelation, meaning that adjacent locations tend to be spatially dispersed (W​o​n​g​ ​&​ ​L​e​e​,​ ​2​0​0​5).

2.4.4 Stationarity

a) Temporal stationarity

Stationarity in a time series means that its statistical properties, such as mean and variance, remain constant over time. A stationarity test was conducted using the Augmented Dickey-Fuller (ADF) test, with hypotheses:

H0: The data is not stationary.

H1: The data is stationary.

ADF is calculated as:

$\Delta Y_t=\beta+\delta Y_{t-1}+\sum_{i=1}^p \alpha_i \Delta Y_{t-i}+\epsilon_t$
(8)

where, $\Delta Y_t$ is the first difference, $\beta$ is the intercept, $\delta$ is the coefficient of the variable to be tested for stationarity, $\Sigma_{i=1}^p \alpha_i \Delta Y_{t-i}$ is the the sum of the first difference lags to the $p$-th lag to capture autocorrelation, and $\epsilon_t$ is the error at time $t$.

If the test statistic is greater than the critical value, H0 cannot be rejected, indicating non-stationary data. To achieve stationarity, differencing was applied (H​y​n​d​m​a​n​ ​&​ ​A​t​h​a​n​a​s​o​p​o​u​l​o​s​,​ ​2​0​1​8):

$y_t^{\prime}=y_t-y_{t-1}$
(9)

where, $y_t^{\prime}$ is the first-order differencing, $y_t$ is the observation at time $t$, and $y_{t-1}$ is the observation at time $(t-$ 1).

b) Spatial stationarity

Stationarity in spatial analysis is an assumption that allows data to be considered as having a uniform level of variation within a region (O​l​i​v​e​r​ ​&​ ​W​e​b​s​t​e​r​,​ ​2​0​1​5). There are three types of stationarity in spatial data analysis:

  • Strict stationarity: Distribution remains the same for all spatial locations.

  • Second-order stationarity: Mean remains constant, and covariance depends only on spatial lag.

  • Intrinsic stationarity: The expectation of differences between locations remains constant.

2.4.5 Identifying the STAR model order

To determine the order of a spatio-temporal model, the spatio-temporal Autocorrelation Function (STACF) and the spatio-temporal Partial Autocorrelation Function (STPACF) can be examined (P​f​e​i​f​e​r​ ​&​ ​D​e​u​t​s​c​h​,​ ​1​9​8​0). STPACF was used to identify the autoregressive component, while STACF was used to identify the moving average component. Table 1 shows the model identification based on the Autocorrelation Function (ACF) and PACF plots.

Table 1. Model identification based on ACF and PACF plots

Model

STACF

STPACF

STAR

tail off

cut off at lag p, spatial lag at λp

STMA

cut off at lag p, spatial lag at λp

tail off

STARMA

tail off

tail off

2.4.6 Parameters estimation

The parameters of the STAR model can be estimated using the Maximum Likelihood Estimation (MLE) method by utilizing the probability density function of ϵ (error), which follows a multivariate normal distribution with a mean of 0 and a constant variance-covariance structure. This was done by defining the likelihood function and identifying the parameter values that maximize the likelihood.

2.4.7 Model diagnostic checking

To ensure the validity of the model used in this study, it is essential to perform diagnostic tests on the residuals. Residuals in a model should satisfy certain assumptions to confirm that the model is correctly specified and provides reliable predictions. These assumptions include white noise properties and normality distribution.

a) Ljung-Box test for white noise residuals

The Ljung-Box test was used to check whether the residuals in the model exhibit the characteristics of white noise. White noise residuals indicate that there is no autocorrelation present in the residuals, meaning that the model has effectively captured the data structure.

The hypotheses for the Ljung-Box test are as follows:

H0: The residuals exhibit white noise properties (no significant autocorrelation).

H1: The residuals do not exhibit white noise properties (presence of autocorrelation).

The test statistic is given by (L​j​u​n​g​ ​&​ ​B​o​x​,​ ​1​9​7​8):

$Q=n(n+2) \sum_{k=1}^K(n-k)^{-1} \hat{\rho}_k^2$
(10)

where, $n$ is the sample size, $\hat{\rho}_k^2$ is the autocorrelation at lag $k$, and $K$ is the number of lags tested.

The test follows a chi-square (χ²) distribution with K degrees of freedom. If the p-value is greater than the significance level (α), then H0 cannot be rejected, indicating that the residuals satisfy the white noise assumption.

b) Mardia’s test for multivariate normality

To verify whether the residuals follow a multivariate normal distribution, Mardia’s test was employed. This test evaluates the skewness and kurtosis of the residuals to assess their normality (M​a​r​d​i​a​,​ ​1​9​7​0). The test statistic for skewness is given by:

$Z_1=\frac{(p+1)(n+1)(n+3)}{6(n+1)(p+1-6)} b_{1, p}$
(11)

where,

$b_{1, p}=\frac{1}{n^2} \sum_{i=1}^n \sum_{j=1}^n\left[\left(x_i-\bar{x}_j\right)^t S^{-1}\left(x_i-\bar{x}_j\right)\right]^3$
(12)

where, $p$ is the number of variables in the multivariate data, $n$ is the sample size, and $S^{-1}$ is the inverse of the covariance matrix.

The test statistic for kurtosis is given by:

$Z_2=\frac{b_{2, p}-p(p+2)}{\sqrt{\frac{8 p(p+2)}{n}}}$
(13)

where,

$b_{2, p}=\frac{1}{n} \sum_{i=1}^n\left[\left(\mathrm{x}_i-\overline{\mathrm{x}}_j\right)^t \mathrm{~S}^{-1}\left(\mathrm{x}_i-\overline{\mathrm{x}}_j\right)\right]^2$
(14)

The hypotheses for Mardia’s test are:

H0: The residuals follow a multivariate normal distribution.

H1: The residuals do not follow a multivariate normal distribution.

If the p-value is greater than the chosen significance level, H0 cannot be rejected, indicating that the residuals follow a normal distribution.

2.4.8 MAPE

MAPE is a statistical measure used to evaluate the accuracy of a forecasting model. It quantifies the percentage error between the actual and predicted values, providing insight into how well the model performs. MAPE is widely used in time-series forecasting as it offers an intuitive and interpretable measure of prediction accuracy.

The formula for MAPE is expressed as (M​a​k​r​i​d​a​k​i​s​ ​e​t​ ​a​l​.​,​ ​1​9​9​8):

$M A P E=\frac{\sum_{t=1}^n\left|P E_t\right|}{n}$
(15)
$P E_t=\frac{y_t-f_t}{y_t} \times 100 \%$
(16)

where, $M A P E$ is the mean absolute percentage error value, $P E_t$ is the percentage error value, $y_t$ is the actual value in period $t, f_t$ is the forecast value in period $t$, and $n$ is the number of observations.

2.4.9 Control chart

In predicting whether a restlessness value is classified as unprecedented or not, this can be determined from the forecast restlessness value that exceeds the UCL on the control chart. The control chart helps in identifying process variability and ensuring that the process remains within the specified control limits.

In this analysis, the focus is only on the UCL because restlessness values exceeding the UCL are the most relevant in identifying a rice crisis. The Lower Control Limit (LCL) is not a concern because values below the LCL are not considered dangerous or indicative of a risk to rice conditions. The calculation of control limits on the control chart is as follows:

$U C L=\bar{y}+L \sigma$
(17)
$C L=\bar{y}$
(18)
$L C L=\bar{y}-L \sigma$
(19)

where, $U C L$ is the upper control limit, $C L$ is the control limit, $L C L$ is the lower control limit, $\bar{y}$ is the mean of all observations, $L$ is the distance for control limits, and $\sigma$ is the standard deviation.

3. Results

3.1 Missing Data Imputation

In this study, there were several missing rice price data for some countries. Missing data can affect analysis results, especially in forecasting models that rely on complete data across locations. Therefore, in this section, the imputation process was conducted to fill in the missing data.

Restlessness values were formed from rice price data, income per capita, and rice consumption. There were no missing data for income per capita and rice consumption, excluding rice price data.

Figure 6. Rice prices (1992-2023)
Figure 7. Rice prices (excluding Brunei)

Figure 6 presents the rice price plot for seven countries (Indonesia, Brunei, Cambodia, Malaysia, the Philippines, Thailand, and Vietnam). It can be seen that Brunei (red line) has a higher price than the other six countries. Brunei's higher rice price makes the plots for the other countries less viable. Therefore, Figure 7 presents the rice price plot for six countries (excluding Brunei).

Based on the rice price plot for six Southeast Asian countries (Indonesia, Cambodia, Malaysia, the Philippines, Thailand, and Vietnam), it is evident that the rice price movements in these six countries have relatively similar patterns, except for Indonesia, where rice prices fluctuated differently from 2011 to 2016 compared to other countries. Overall, the price trends tend to move in the same direction with closely related fluctuations among countries, indicating a common spatial pattern.

Therefore, this study performs imputation using spatial interpolation with the IDW method, following the formula explained in subsection 2.2. This method assigns weights based on the distance between locations, where closer locations have a greater influence on estimating the missing values.

Figure 8. Rice prices for seven countries after imputation

Figure 8 represents the rice price plot for seven countries (Indonesia, Brunei, Cambodia, Malaysia, the Philippines, Thailand, and Vietnam) after imputation using IDW. Figure 9 shows the rice price plot for six countries (Indonesia, Cambodia, Malaysia, the Philippines, Thailand, and Vietnam) after imputation. The existing data for Brunei's rice prices consistently show significantly higher values than the other six countries.

Figure 9. Rice prices for six countries (excluding Brunei) after imputation

After spatial imputation, Brunei's imputed rice prices were distorted downward, making the pattern more similar to that of other countries. The distortion of Brunei's rice price imputation may reduce the validity of the analysis results.

Figure 10. Combined rice price control chart

Figure 10 displays the control chart for the combined rice prices of each country. Brunei has rice prices that exceed the UCL and significantly differ from those of other countries. This indicates that Brunei's rice prices do not resemble those of other countries and, therefore, cannot be forced into spatial imputation. Consequently, Brunei was excluded from the main analysis. Special handling is required to impute Brunei's rice prices to ensure more accurate results that align with the actual data characteristics.

3.2 Descriptive Statistics

In this study, the data used consists of annual restlessness values from the period of 1992–2023 in Southeast Asia, covering six countries: Indonesia, the Philippines, Thailand, Vietnam, Malaysia, and Cambodia. The summary of central tendency and dispersion measures of the data is presented in Table 2.

Table 2. Descriptive statistics of restlessness in each country

Country

Min

Max

Median

Mean

Variance

St. Dev.

Indonesia

0.0094

0.0397

0.025865

0.0247

0.0000875

0.0093588

Cambodia

0.0461

0.13542

0.0691

0.0718

0.0004275

0.0206764

Malaysia

0.0019

0.004023

0.0024

0.0027

0.0000004

0.0006470

Philippines

0.0125

0.02630

0.0192

0.0188

0.0000126

0.0035565

Thailand

0.0059

0.01243

0.0083

0.008467

0.0000025

0.0015901

Vietnam

0.0164

0.16506

0.0514

0.0595

0.0018961

0.0435444

Based on Table 2, the country with the highest average restlessness value is Vietnam at 0.0595, while the lowest is Malaysia at 0.0027. This indicates that Vietnam has experienced greater fluctuations in rice-related economic and social conditions compared to Malaysia. Furthermore, Cambodia has the highest variance (0.0004275), showing that its restlessness values have fluctuated significantly over time. In contrast, Malaysia has the lowest variance (0.0000004), indicating a relatively stable trend in restlessness values. The standard deviation values confirm these findings, with Vietnam and Cambodia having the highest standard deviations (0.0435444 and 0.0206764, respectively), reflecting higher unpredictability in their restlessness indices. Conversely, Malaysia and Thailand have the lowest standard deviations (0.0006470 and 0.0015901, respectively), indicating more stable trends in restlessness over the years.

To further examine the restlessness phenomenon over time, it is important to analyze how these trends evolve temporally within each country. Thus, the following plots illustrate the time series trends of restlessness values across different countries.

Figure 11. Restlessness plot

Figure 11 presents the restlessness plot for six Southeast Asian countries: Indonesia, Cambodia, Malaysia, the Philippines, Thailand, and Vietnam. It reveals a general downward trend in restlessness levels from 1992 to 2023. Vietnam and Cambodia initially exhibited relatively high restlessness levels in the 1990s, but these gradually declined over time.

As shown in Figure 12, the Philippines, Indonesia, and Thailand demonstrate relatively stable restlessness levels throughout the observed period. Malaysia, on the other hand, displays the most stable and lowest restlessness values with minimal fluctuations. Overall, these trends indicate that, despite fluctuations in certain periods, most countries exhibit a declining pattern in restlessness levels over time.

Figure 12. Restlessness plot (excluding Cambodia and Vietnam)
3.3 Spatial Weighting

The values of the inverse distance weights were obtained from calculations based on the actual distances between locations, considering their latitude and longitude coordinates. The latitude and longitude coordinates of each country are presented in Table 3.

Table 3. Latitude and longitude coordinates

Country

Latitude

Longitude

Indonesia

-0.789275

113.921327

Cambodia

12.565679

104.990963

Malaysia

4.210484

101.975766

Philippines

12.879721

121.774017

Thailand

15.870032

100.992541

Vietnam

14.058324

108.277199

The latitude and longitude coordinates were then converted into kilometers so that the resulting inverse distance weight matrix calculations, using Eq. (5), are displayed in the following matrix:

$W=\left[\begin{array}{cccccc}0 & 0.199 & 0.246 & 0.202 & 0.152 & 0.201 \\ 0.087 & 0 & 0.157 & 0.085 & 0.274 & 0.396 \\ 0.189 & 0.276 & 0 & 0.115 & 0.209 & 0.210 \\ 0.214 & 0.206 & 0.158 & 0 & 0.166 & 0.256 \\ 0.092 & 0.380 & 0.166 & 0.095 & 0 & 0.267 \\ 0.098 & 0.439 & 0.133 & 0.118 & 0.213 & 0\end{array}\right]$

3.4 Spatial Autocorrelation Testing (MoranST)

At this stage, spatial autocorrelation testing was conducted using MoranST. The output from MoranST with the assistance of R software is shown in Table 4.

Table 4. MoranST spatial autocorrelation test results

MoranST Value

p-value

0.6000271

0.00990099

The test results shown in Table 4 include a MoranST value of 0.6000271 and a p-value of 0.00990099. Using a significance level (α) of 5%, since the p-value < 0.05, the decision taken is to reject H0 in favor of the alternative hypothesis H1. This means that there is significant spatial autocorrelation. Practically, this indicates that neighboring regions tend to have similar characteristics.

3.5 Data Training and Testing Split

In this study, the restlessness data were divided into two parts: in-sample and out-of-sample data. The training data, or in-sample data, was used for model building. Meanwhile, the testing data, or out-of-sample data, was used to evaluate the model’s predictive capability.

Various splits between training and testing data (in-sample and out-of-sample) were tested to assess the model’s performance in generating optimal forecasts. These variations aim to identify the combination that provides the best results based on the predefined evaluation criteria, in this case, using MAPE as the primary indicator.

However, due to the large number of tested variations, not all results can be presented in detail in this chapter. Therefore, the analysis in this chapter focuses on the data split variation that yields the best MAPE value, specifically the 80% in-sample and 20% out-of-sample split. The 80% in-sample data cover the period from 1992 to 2017 for model training and development. Meanwhile, the 20% out-of-sample data covers the period from 2018 to 2023 for testing the model’s predictive capability on unseen data.

Figure 13. In-sample and out-of-sample data division for all analyzed countries

Figure 13 shows the plot of the results of dividing the in-sample and out-of-sample data for the six countries analyzed. The solid line represents the in-sample data, which is used to train and build the model, while the dashed line represents the out-of-sample data, which is used to test the predictive ability of the model. Different colors of the lines indicate different country representations. Figure 14 shows the plot of the results of dividing the in-sample and out-of-sample data for each country separately.

Figure 14. In-sample and out-of-sample data division for each country separately
3.6 Stationarity Test

The ADF test results for the training data on restlessness in six countries are as follows:

Table 5. ADF training data for restlessness in six countries

Country

ADF Test Before Differencing

ADF Test After 1st Differencing

ADF Test After 2nd Differencing

ADF Test After 3rd Differencing

p-value

p-value

p-value

p-value

Indonesia

0.5297 (non-stationarity)

0.3403 (non-stationarity)

0.01 (stationary)

0.01 (stationary)

Cambodia

0.2032 (non-stationarity)

0.01098 (stationary)

0.0916 (non-stationarity)

0.01 (stationary)

Malaysia

0.4488 (non-stationarity)

0.01 (stationary)

0.01 (stationary)

0.01 (stationary)

Philippines

0.2999 (non-stationarity)

0.1109 (non-stationary)

0.0476 (stationary)

0.01 (stationary)

Thailand

0.2755 (non-stationarity)

0.01641 (stationary)

0.01 (stationary)

0.01 (stationary)

Vietnam

0.7064 (non-stationarity)

0.2914 (non-stationarity)

0.0794 (non-stationarity)

0.02372 (stationary)

Based on Table 5, it can be seen that before differencing, the restlessness data for all six countries had a p-value greater than 0.05. If the p-value < 0.05, the significance level, then the null hypothesis is rejected, concluding that the data is stationary. Therefore, the restlessness data of the six countries was non-stationary because the p-value was greater than 0.05, and differencing was necessary to remove trends or seasonal patterns causing the non-stationarity. After performing differencing three times, a p-value < 0.05 was obtained, indicating that the data had become stationary.

3.7 Identification of the STAR Model Order

The spatial order used in this study is spatial order 1 because each research location is within the same region, namely Southeast Asia. The identification of the autoregressive and moving average temporal orders was carried out by examining the STACF and STPACF plots generated from differenced data using IDW.

Figure 15. STACF plot
Figure 16. STPACF plot

Based on Figure 15, the STACF plot shows no significant lag in spatial lag = 1. Therefore, there is no moving average component in this modeling. Furthermore, based on Figure 16, the STPACF plot indicates significant lags at lags 4, 5, and 8. Hence, the suitable candidate models are STAR models with AR orders of 4, 5, or 8.

3.7.1 Selecting the best model

Next, the best model was selected from several candidate models by examining the Bayesian Information Criterion (BIC) value for each model. BIC is a statistical measure used for model selection from a set of candidate models. The BIC value is based on the likelihood value and a complexity penalty. A lower BIC value indicates a lower penalty, making the model preferable. The BIC values for each candidate model are presented in Table 6.

Table 6. BIC of candidate models

Model

BIC

AR = 4

-796.1193

AR = 5

-791.1394

AR = 8

-678.3848

AR indicates autoregressive.

Table 6 shows the BIC values for each candidate model. The lowest BIC value is found in the AR = 4 model with a BIC of -796.1193. However, since the difference between the BIC values of AR = 4 and AR = 5 is small, the next step is to compare the Root Mean Squared Error (RMSE) and Mean Squared Error (MSE) of the residuals for these two models.

Table 7. RMSE and MSE of candidate models

Model

RMSE

MSE

AR = 4

0.0165414

0.0002736179

AR = 5

0.01629915

0.0002656623

Note: AR indicates autoregressive.

As shown in Table 7, it can be observed that the model with AR order 5 has lower RMSE and MSE values compared to AR order 4. Therefore, the selected model is STAR (5,1).

3.8 Parameter Estimation

Parameter estimation was carried out based on the location weights used. The estimation method employed is MLE. The calculation results, with the assistance of R software, yield the estimated parameter values as follows:

$\begin{aligned} \hat{\phi}_{10} & =0.914009, \hat{\phi}_{11}=0.1508, \hat{\phi}_{20}=0.051005, \hat{\phi}_{21}=-0.0839, \hat{\phi}_{30}=0.117822, \hat{\phi}_{31}= \\ & -0.1703, \hat{\phi}_{40}=-0.223924, \hat{\phi}_{41}=0.1057, \hat{\phi}_{50}=0.073784, \operatorname{and} \hat{\phi}_{51}=0.0163 .\end{aligned}$

The parameters of the STAR (5,1) model indicate the magnitude of the influence of a location on itself and the effect of a location on its surrounding locations from lag 1 to lag 5. The estimated STAR (5,1) model, constructed using the inverse distance weight matrix, is as follows:

$\left[\begin{array}{l}\widehat{Z_1}(t) \\ \widehat{Z_2}(t) \\ \widehat{Z_3}(t) \\ \widehat{Z_4}(t) \\ \widehat{Z_5}(t) \\ \widehat{Z_6}(t)\end{array}\right]=0.914009\left[\begin{array}{l}\widehat{Z_1}(t-1) \\ \widehat{Z_2}(t-1) \\ \widehat{Z_3}(t-1) \\ \widehat{Z_4}(t-1) \\ \widehat{Z_5}(t-1) \\ \widehat{Z_6}(t-1)\end{array}\right]$$+0.1508\left[\begin{array}{cccccc}0 & 0.199 & 0.246 & 0.202 & 0.152 & 0.201 \\ 0.087 & 0 & 0.157 & 0.085 & 0.274 & 0.396 \\ 0.189 & 0.276 & 0 & 0.115 & 0.209 & 0.210 \\ 0.214 & 0.206 & 0.158 & 0 & 0.166 & 0.256 \\ 0.092 & 0.380 & 0.166 & 0.095 & 0 & 0.267 \\ 0.098 & 0.439 & 0.133 & 0.118 & 0.213 & 0\end{array}\right]\left[\begin{array}{l}\widehat{Z}_1(t-1) \\ \widehat{Z}_2(t-1) \\ \widehat{Z}_3(t-1) \\ \widehat{Z}_4(t-1) \\ \widehat{Z}_5(t-1) \\ \widehat{Z}_6(t-1)\end{array}\right]+0.051005\left[\begin{array}{l}\widehat{Z}_1(t-2) \\ \widehat{Z_2}(t-2) \\ \widehat{\widehat{Z}_3}(t-2) \\ \widehat{Z_4}(t-2) \\ \widehat{Z}_5(t-2) \\ \widehat{Z}_6(t-2)\end{array}\right]$

$-0.0839\left[\begin{array}{cccccc}0 & 0.199 & 0.246 & 0.202 & 0.152 & 0.201 \\ 0.087 & 0 & 0.157 & 0.085 & 0.274 & 0.396 \\ 0.189 & 0.276 & 0 & 0.115 & 0.209 & 0.210 \\ 0.214 & 0.206 & 0.158 & 0 & 0.166 & 0.256 \\ 0.092 & 0.380 & 0.166 & 0.095 & 0 & 0.267 \\ 0.098 & 0.439 & 0.133 & 0.118 & 0.213 & 0\end{array}\right]\left[\begin{array}{c}\widehat{Z}_1(t-2) \\ \widehat{Z}_2(t-2) \\ \widehat{Z}_3(t-2) \\ \widehat{Z}_4(t-2) \\ \widehat{Z}_5(t-2) \\ \widehat{Z}_6(t-2)\end{array}\right]+0.117822\left[\begin{array}{l}\widehat{Z}_1(t-3) \\ \widehat{Z}_2(t-3) \\ \widehat{Z}_3(t-3) \\ \widehat{Z}_4(t-3) \\ \widehat{Z_5}(t-3) \\ \widehat{Z_6}(t-3)\end{array}\right]-$

$0.1703\left[\begin{array}{cccccc}0 & 0.199 & 0.246 & 0.202 & 0.152 & 0.201 \\ 0.087 & 0 & 0.157 & 0.085 & 0.274 & 0.396 \\ 0.189 & 0.276 & 0 & 0.115 & 0.209 & 0.210 \\ 0.214 & 0.206 & 0.158 & 0 & 0.166 & 0.256 \\ 0.092 & 0.380 & 0.166 & 0.095 & 0 & 0.267 \\ 0.098 & 0.439 & 0.133 & 0.118 & 0.213 & 0\end{array}\right]\left[\begin{array}{l}\widehat{Z_1}(t-3) \\ \widehat{Z_2}(t-3) \\ \widehat{Z_3}(t-3) \\ \widehat{Z_4}(t-3) \\ \widehat{Z_5}(t-3) \\ \widehat{Z_6}(t-3)\end{array}\right]-0.223924\left[\begin{array}{l}\widehat{Z_1}(t-4) \\ \widehat{Z_2}(t-4) \\ \widehat{Z_3}(t-4) \\ \widehat{Z_4}(t-4) \\ \widehat{Z_5}(t-4) \\ \widehat{Z_6}(t-4)\end{array}\right]+$

$0.1057\left[\begin{array}{cccccc}0 & 0.199 & 0.246 & 0.202 & 0.152 & 0.201 \\ 0.087 & 0 & 0.157 & 0.085 & 0.274 & 0.396 \\ 0.189 & 0.276 & 0 & 0.115 & 0.209 & 0.210 \\ 0.214 & 0.206 & 0.158 & 0 & 0.166 & 0.256 \\ 0.092 & 0.380 & 0.166 & 0.095 & 0 & 0.267 \\ 0.098 & 0.439 & 0.133 & 0.118 & 0.213 & 0\end{array}\right]\left[\begin{array}{l}\widehat{Z_1}(t-4) \\ \widehat{Z_2}(t-4) \\ \widehat{Z_3}(t-4) \\ \widehat{Z_4}(t-4) \\ \widehat{Z_5}(t-4) \\ \widehat{Z_6}(t-4)\end{array}\right]+0.073784\left[\begin{array}{l}\widehat{Z_1}(t-5) \\ \widehat{Z_2}(t-5) \\ \widehat{Z_3}(t-5) \\ \widehat{Z_4}(t-5) \\ \widehat{Z_5}(t-5) \\ \widehat{Z_6}(t-5)\end{array}\right]+$

$0.0163\left[\begin{array}{cccccc}0 & 0.199 & 0.246 & 0.202 & 0.152 & 0.201 \\ 0.087 & 0 & 0.157 & 0.085 & 0.274 & 0.396 \\ 0.189 & 0.276 & 0 & 0.115 & 0.209 & 0.210 \\ 0.214 & 0.206 & 0.158 & 0 & 0.166 & 0.256 \\ 0.092 & 0.380 & 0.166 & 0.095 & 0 & 0.267 \\ 0.098 & 0.439 & 0.133 & 0.118 & 0.213 & 0\end{array}\right]\left[\begin{array}{l}\widehat{Z}_1(t-5) \\ \widehat{Z}_2(t-5) \\ \widehat{Z}_3(t-5) \\ \widehat{Z}_2(t-5) \\ \widehat{Z}_5(t-5) \\ \widehat{Z}_6(t-5)\end{array}\right]$

3.9 Model Diagnostic Check

The model evaluation was conducted by examining the white noise assumption and normality of the residuals.

Table 8 presents the results of the Ljung-Box test. The test results indicate that all p-values are greater than 0.05 for each lag, meaning that the model fails to reject the null hypothesis (H0). In other words, there is no autocorrelation among the residuals, confirming that the residuals meet the white noise assumption.

As shown in Table 9, the p-values for both skewness and kurtosis are greater than 0.05, leading to the conclusion that the residuals follow a multivariate normal distribution, and the model is valid for use.

Table 8. Ljung-Box test

Lags

Statistic

df

p-value

1

43.49246

36

0.1825634

2

75.22718

72

0.3743311

3

105.02482

108

0.5631022

4

132.94437

144

0.7354454

5

152.27062

180

0.9342813

6

178.37923

216

0.9709345

7

199.95495

252

0.9932467

8

223.65669

288

0.9980316

9

244.10823

324

0.9996806

10

265.99832

360

0.9999384

11

282.3011

396

0.9999963

12

299.51967

432

0.9999998

13

314.79817

468

1

14

328.87819

504

1

15

339.19039

540

1

16

350.80068

576

1

17

360.29296

612

1

18

364.58606

648

1

19

367.9071

684

1

20

369.96748

720

1

Table 9. Multivariate normality test

Test

Statistic

p-value

Result

Mardia skewness

66.82609927

0.152515769

Yes

Mardia kurtosis

0.1601768381

0.8727417743

Yes

3.10 MAPE

Out-of-sample forecasting is useful for evaluating the model’s ability to predict future data. The model’s accuracy was assessed using the MAPE metric. Table 10 presents the out-of-sample forecasting results of restlessness values in the six countries.

It can be observed that the MAPE values for the out-of-sample restlessness data for each country vary. The smaller the MAPE value, the more accurate the forecast. The overall MAPE value for all countries is 15.1%. Vietnam, Indonesia, the Philippines, and Thailand have excellent MAPE values, all below 15%. Meanwhile, Cambodia and Malaysia have slightly higher MAPE values, at 19% and 36%, respectively. Below are the forecast restlessness values using the STAR (5,1) model presented in a time series plot.

Table 10. Out-of-sample forecasting

Country

Year

Actual Value

Forecast Value

MAPE

Indonesia

2018

0.013690239

0.01270866

7.167575%

2019

0.012292445

0.01193459

2020

0.009925645

0.01189298

2021

0.010678892

0.00995473

2022

0.009997526

0.00943791

2023

0.009355219

0.00959569

Cambodia

2018

0.0492267

0.05408268

19.03478%

2019

0.04540735

0.05061554

2020

0.0419276

0.05119968

2021

0.03925779

0.05267709

2022

0.03658789

0.04615783

2023

0.03403994

0.04901998

Malaysia

2018

0.002119712

0.00227947

36.11504%

2019

0.002281153

0.00230061

2020

0.002953866

0.0024535

2021

0.003152864

0.00224657

2022

0.003411541

0.00195969

2023

0.003641943

0.00209289

Philippines

2018

0.01301014

0.01559437

10.62152%

2019

0.01250111

0.01520368

2020

0.01244133

0.01381602

2021

0.01178882

0.01329812

2022

0.01135831

0.01274531

2023

0.0109332

0.01276547

Thailand

2018

0.006160782

0.00713233

13.70015%

2019

0.006271137

0.00774868

2020

0.006580454

0.00728447

2021

0.006569694

0.00793326

2022

0.006702965

0.00760609

2023

0.00668659

0.00759141

Vietnam

2018

0.02035355

0.01967505

4.04275%

2019

0.01919096

0.0197425

2020

0.01833697

0.01829328

2021

0.01759828

0.01969642

2022

0.01675452

0.01662018

2023

0.01598074

0.01640847

Overall MAPE

15.11364

Figure 17 illustrates the forecast values generated by the STAR (5,1) model for the out-of-sample data. The model’s predictions (blue line) generally follow the pattern of the actual values (red line) fairly well, except for Malaysia. The forecast results for Malaysia exhibit higher errors compared to other countries, as indicated by its MAPE, which is significantly higher at 36% compared to other countries. The out-of-sample forecasting for Malaysia shows a trend of increasing restlessness values, whereas the actual values do not exhibit this trend.

Additionally, Malaysia's restlessness values are much lower, ranging between 0.001 and 0.004, whereas other countries have restlessness values ranging from 0.01 to 0.1. The STAR model that was constructed pulls Malaysia’s forecasts upward due to the higher restlessness values in other countries, as the model assumes all locations to be homogeneous.

Table 11 is a summary of experiments using different data splits and model orders.

Figure 17. Actual vs. out-of-sample forecast data
Table 11. Summary of data splitting variations and models

No.

Data Split

Model

Diagnostic Check

MAPE

In-sample

Out-of-sample

1

60/40

STAR (10,1)

Residuals are not normal.

-

-

2

75/25

STAR (5,1)

Residuals are not normal.

-

-

3

75/25

STAR (4,1)

Valid

15%

31%

4

70/30

STAR (3,1)

Valid

16%

46%

5

80/20

STAR (5,1)

Valid

15.69%

15.11%

6

80/20

STAR (3,1)

Valid

15.60%

19.30%

7

80/20

STAR (4,1)

Valid

15%

19.19%

8

90/10

STAR (6,1)

Valid

17%

26.22%

3.11 Forecasting

Next, the STAR (5,1) model was used to forecast restlessness values over the next five years, from 2024 to 2028. This forecasting process utilizes all available data without separating it into in-sample and out-of-sample data. Figure 18 presents the five-year forecasting results of restlessness values.

Figure 18. Five-year forecast of restlessness

Figure 18 illustrates the restlessness forecasts for the six analyzed countries over the next five years. The model’s predictions provide insights into potential future trends in restlessness values, which may be useful for policy-making and crisis-prevention strategies.

3.12 Control Chart

Next, to determine whether a country in a given year falls into an unprecedented condition or not, a further analysis was conducted using a control chart. The control chart includes two control limits, namely the UCL and the LCL. The UCL line serves as a threshold to identify the potential rice crisis in the country. If the restlessness value exceeds the UCL line, the country is categorized as experiencing unprecedented restlessness, indicating a high likelihood of a rice crisis. Conversely, if the value remains below the UCL, the restlessness condition is considered to be within an acceptable control limit.

The control chart calculation utilizes the z-score distance for the UCL and LCL. In this study, the control limits were established when the data deviated by 90% from the mean. Based on the z-score table, the control limit distance (L) for 90% is 1.29. Thus, the calculation of UCL and LCL is as follows:

$U C L=\bar{y}+1.29 \sigma$
(20)
$U C L=\bar{y}-1.29 \sigma$
(21)
Figure 19. Control chart of restlessness values

Figure 19 illustrates the control chart for six countries, namely Indonesia, Cambodia, Malaysia, the Philippines, Thailand, and Vietnam. In each plot, the blue line represents the actual restlessness values, while the red line represents the predicted restlessness values for the next five years.

Based on the figure, the forecast values for the next five years in all countries remain below the UCL. Thus, no country is detected as experiencing unprecedented restlessness or a rice crisis in the coming five years.

4. Discussion

The focus of this research stems from concerns about the rice crisis in Southeast Asia due to various factors such as climate change, strict export policies, and fluctuations in rice prices. One significant event that triggered this study was India's rice export ban in 2023, which caused a surge in rice prices in several countries. With the potential for a future rice crisis, a method is needed to forecast possible crises that may occur. Therefore, this study uses the restlessness indicator, calculated based on rice prices, rice consumption, and per capita income, to assess the potential for food instability. The STAR model was used to forecast restlessness values in six Southeast Asian countries, namely Indonesia, the Philippines, Thailand, Vietnam, Malaysia, and Cambodia, using secondary data from 1992 to 2023.

Among the various candidate models tested, the STAR (5,1) model was selected as the best model based on model selection criteria such as the BIC, RMSE, and MSE. This model considers temporal dependencies up to five years back and spatial dependencies up to one level of inter-country relationships. The research data was divided into training data (1992-2018) and testing data (2019-2023) to evaluate the model's performance in forecasting restlessness values.

To ensure the validity of the model, diagnostic tests were conducted, including the white noise test, which showed that the residuals of the model met the white noise assumption without systematic patterns in the model's residual predictions. The normality test using Mardia's test indicated that the residuals followed a multivariate normal distribution. Additionally, Moran’s I test confirmed that there was a significant spatial relationship in the data, supporting the selection of the STAR model.

The control chart was used to identify whether unprecedented restlessness occurred, serving as an indicator of a potential rice crisis in the next five years (2024-2028). The analysis results showed that no country experienced unprecedented restlessness during the forecast period. The restlessness values over the next five years remained within control limits, indicating a low potential for a rice crisis. Although there was a slight increase in restlessness in some countries, the values remained within a reasonable range and did not show signs of extreme instability.

Based on the results of this study, several measures can be taken, including regular monitoring of restlessness. Although there is no indication of a crisis in the next five years, restlessness values should continue to be observed to anticipate any significant increases early. Food source diversification is also necessary for Southeast Asian countries to develop alternative food sources besides rice to reduce dependency on a single commodity. Furthermore, regional cooperation in rice trade and production can help maintain price and supply stability. Improving production efficiency through investments in irrigation technology and precision agriculture can also help secure future rice supplies.

For the Indonesian government, several aspects of cooperation with the Association of Southeast Asian Nations (ASEAN) countries in the rice sector should be considered:

First, expanding the export market and diversifying rice-based products with Vietnam. Indonesia and Vietnam can collaborate to create high-value-added rice-based products such as rice flour, cereals, or processed foods. This initiative aims to expand ASEAN and global export markets while enhancing the competitiveness of both countries.

Second, improving rice logistics efficiency and food diversification with the Philippines. The Philippines can work with Indonesia to establish direct trade routes to lower distribution costs. In non-crisis conditions, both countries can design flexible import and export policies to prevent overstocking and optimize surplus production.

Third, developing agricultural technology and ASEAN rice safety standards with Cambodia and Malaysia. Under stable conditions, Indonesia can assist Cambodia and Malaysia in developing digital agricultural technologies such as the Internet of Things (IoT), drone monitoring, and smart irrigation to improve production efficiency. Additionally, cooperation in establishing ASEAN rice quality standards can be strengthened to position ASEAN as a premium rice exporter in the global market.

Finally, stabilizing prices and strengthening rice reserves with Thailand. In the absence of a rice crisis threat, both countries can develop a more adaptive price stabilization system and establish regional rice reserve policies to anticipate economic uncertainties or natural disasters that may disrupt future food security.

This study demonstrates that the STAR (5,1) model is capable of forecasting restlessness values with reasonably good accuracy, as indicated by a MAPE value of 15.1%. The forecast results show no indication of unprecedented restlessness, meaning there are no significant signs of a rice crisis in the next five years. However, the previously mentioned policy measures still need to be implemented to ensure future food stability.

5. Conclusions

Based on the results and discussion above, the best-performing model obtained in this study is the STAR(5,1) model. The general form of the STAR model is expressed as:

$y(t)=\sum_{k=1}^5 \Sigma_{l=0}^1\left[\phi_{k l} W^{(l)} y(t-k)\right]+\epsilon(t)$
(22)

where, $W^{(l)}$ represents the spatial weight matrix, $\phi_{k l}$ is the estimated parameters, and $\epsilon$ is the error term. The spatial weight matrix used in this study is based on inverse distance, which accounts for the geographical proximity between countries. The estimated parameters indicate that both temporal and spatial dependencies contribute to the variations in restlessness values across Southeast Asian nations.

The results of the control chart analysis show that no country is classified as experiencing unprecedented restlessness in the five-year forecast period from 2024 to 2028. This suggests that, under current economic and agricultural conditions, the likelihood of a severe rice crisis remains low. While minor fluctuations in restlessness values were observed, they remained within control limits, indicating stable rice market conditions.

The recommendation for future research is to develop the model by adding a dummy variable that represents the difference in the range of restlessness values for Malaysia. The difference in the range of restlessness values in Malaysia compared to other countries may affect the model's performance in forecasting. By incorporating a dummy variable, the model can be more adaptive in capturing the differences in the range of restlessness values occurring in Malaysia, thereby improving the accuracy of predictions for the country. Additionally, further development can be carried out by using consumption data without missing values and by seeking consumption data at the consumer level.

Based on the findings of this research, several policy recommendations can be proposed. Diversification of rice supply sources should be encouraged to reduce reliance on specific suppliers. Strengthening regional cooperation in rice trade and production can enhance supply chain resilience, and investments in agricultural technology and irrigation infrastructure can improve productivity and sustainability. By implementing these strategies, Southeast Asian nations can mitigate potential risks and ensure stability in the rice market.

Author Contributions

Conceptualization, Y.H. and C.D.A.M.; Methodology, C.D.A.M., Y.H. and B.T.; software, C.D.A.M.; validation, Y.H. and B.T.; formal analysis, Y.H. and B.T.; investigation, C.D.A.M. and Y.H.; resources, Y.H. and B.T.; data curation, C.D.A.M.; writing—original draft preparation, C.D.A.M.; writingreview and editing, C.D.A.M., Y.H. and B.T.; visualization, C.D.A.M.; supervision, Y.H. and B.T.; project administration, Y.H.; funding acquisition, Y.H. All authors have read and agreed to the published version of the manuscript.

Funding
The paper was funded by the Directorate of Research and Community Engagement of Padjadjaran University.
Data Availability

The data (rice price, income per capita, and rice consumption) supporting our research results are openly available and can be accessed from FAO (https://www.fao.org/faostat/en/#data/PP), World Bank (https://data.worldbank.org/indicator/NY.GDP.PCAP.CD), and FAS (https://apps.fas.usda.gov/psdonline/app/index.html#/app/advQuery).

Acknowledgments

The authors would like to express their sincere gratitude to the Directorate of Research and Community Engagement of Padjadjaran University for their valuable support during the writing of this paper and for providing financial support for publication in this journal.

Conflicts of Interest

The authors declare no conflict of interest.

References
Anderson, C. & Ryan, L. (2017). A comparison of spatio-temporal disease mapping approaches including an application to ischaemic heart disease in New South Wales, Australia. Int. J. Environ. Res. Public Health, 14(2), 146. [Google Scholar] [Crossref]
FAO. (2023). The State of Food Security and Nutrition in the World 2023. Rome, Italy: FAO; IFAD; UNICEF; WFP; WHO. [Google Scholar] [Crossref]
Hidayat, Y., Purwandari, T., Ratnasari, D., Sukono, Saputra, J., & Subiyanto. (2021). Objective method for determining the importance of unprecedented restlessness as a rice crisis indicator at the national level. Agronomy, 11(6), 1195. [Google Scholar] [Crossref]
Hyndman, R. J. & Athanasopoulos, G. (2018). Forecasting: Principles and Practice (2nd edition). OTexts: Melbourne, Australia. https://otexts.com/fpp2/ [Google Scholar]
International Research Institute for Climate & Society. (2016). Map of El Nino. https://iridl.ldeo.columbia.edu/maproom/IFRC/FIC/elninorain.html#tabs-2 [Google Scholar]
Kano, K., Kano, T., & Takechi, K. (2022). The price of distance: Pricing-to-market and geographic barriers. J. Econ. Geogr., 22(4), 873–899. [Google Scholar] [Crossref]
Kim, S. (2024). Transmission of economic growth between neighboring countries in Europe. Appl. Econ., 1–16. [Google Scholar] [Crossref]
LeSage, J. P. (2008). An introduction to spatial econometrics. Rev. Econ. Ind., 19–44. [Google Scholar] [Crossref]
Ljung, G. M. & Box, G. E. P. (1978). On a measure of lack of fit in time series models. Biometrika, 65(2), 297–303. [Google Scholar] [Crossref]
Lu, G. Y. & Wong, D. W. (2008). An adaptive inverse-distance weighting spatial interpolation technique. Comput. Geosci., 34(9), 1044–1055. [Google Scholar] [Crossref]
Makridakis, S., Wheelwright, S. C., & Hyndman, R. J. (1998). Forecasting: Methods and Applications (Third Edition). New York: John Wiley & Sons. [Google Scholar]
Mardia, K. V. (1970). Measures of multivariate skewness and kurtosis with applications. Biometrika, 57(3), 519–530. [Google Scholar] [Crossref]
Mustafa, M., Yuyun, H., & Ismail, M. (2015). Developing unprecedented restlessness as the new strong indicator of rice crisis at national level. Glob. J. Pure Appl. Math., 11(6), 4139–4160. [Google Scholar]
Oliver, M. A. & Webster, R. (2015). Basic Steps in Geostatistics: The Variogram and Kriging. Springer. [Google Scholar] [Crossref]
Pfeifer, P. E. & Deutsch, S. J. (1980). A three-stage iterative procedure for space-time modeling. Technometrics, 22(1), 35–47. [Google Scholar] [Crossref]
Rahmawati, S. (2006). Status perkembangan perbaikan sifat genetik padi menggunakan transformasi argobacterium. J. Agrobiogen, 2(1), 36–44. [Google Scholar]
Ruchjana, B. N. (2019). Pengembangan model spatio temporal dan aplikasinya. In Prosiding Seminar Nasional Matematika Dan Pendidikan Matematika, UIN Raden Intan Lampung, 2(1). https://ejournal.radenintan.ac.id/index.php/pspm/article/view/4295 [Google Scholar]
Scaife, A., Guilyardi, E., Cain, M., & Gilbert, A. (2019). What is the El Niño–Southern Oscillation? Weather, 74(7), 250–251. [Google Scholar] [Crossref]
Soni, A. & Singh, P. (2022). A review paper on El Nino. Int. J. Innov. Res. Eng. Manag., 9(1), 269–272. [Google Scholar] [Crossref]
Tobler, W. R. (1970). A computer movie simulating urban growth in the Detroit region. Econ. Geogr., 46, 234–240. [Google Scholar] [Crossref]
Wong, D. W. S. & Lee, J. (2005). Statistical Analysis of Geographic Information with ArcView GIS and ArcGIS. John Wiley & Sons, Inc. [Google Scholar]

Cite this:
APA Style
IEEE Style
BibTex Style
MLA Style
Chicago Style
GB-T-7714-2015
Hidayat, Y., Tantular, B., & Manurung, C. D. A. (2025). Forecasting Rice Crisis Dynamics in Southeast Asia Using a Spatio-Temporal Autoregressive Model Based on the Restlessness Indicator. Org. Farming, 11(1), 13-38. https://doi.org/10.56578/of110102
Y. Hidayat, B. Tantular, and C. D. A. Manurung, "Forecasting Rice Crisis Dynamics in Southeast Asia Using a Spatio-Temporal Autoregressive Model Based on the Restlessness Indicator," Org. Farming, vol. 11, no. 1, pp. 13-38, 2025. https://doi.org/10.56578/of110102
@research-article{Hidayat2025ForecastingRC,
title={Forecasting Rice Crisis Dynamics in Southeast Asia Using a Spatio-Temporal Autoregressive Model Based on the Restlessness Indicator},
author={Yuyun Hidayat and Bertho Tantular and Cosmas David Ananda Manurung},
journal={Organic Farming},
year={2025},
page={13-38},
doi={https://doi.org/10.56578/of110102}
}
Yuyun Hidayat, et al. "Forecasting Rice Crisis Dynamics in Southeast Asia Using a Spatio-Temporal Autoregressive Model Based on the Restlessness Indicator." Organic Farming, v 11, pp 13-38. doi: https://doi.org/10.56578/of110102
Yuyun Hidayat, Bertho Tantular and Cosmas David Ananda Manurung. "Forecasting Rice Crisis Dynamics in Southeast Asia Using a Spatio-Temporal Autoregressive Model Based on the Restlessness Indicator." Organic Farming, 11, (2025): 13-38. doi: https://doi.org/10.56578/of110102
HIDAYAT Y, TANTULAR B, MANURUNG C D A. Forecasting Rice Crisis Dynamics in Southeast Asia Using a Spatio-Temporal Autoregressive Model Based on the Restlessness Indicator[J]. Organic Farming, 2025, 11(1): 13-38. https://doi.org/10.56578/of110102
cc
©2025 by the author(s). Published by Acadlore Publishing Services Limited, Hong Kong. This article is available for free download and can be reused and cited, provided that the original published version is credited, under the CC BY 4.0 license.