20250812 meeting

Introduction

Due to the combined effects of human activities and natural climate variability (e.g., El Niño), global surface temperatures have shown a long-term upward trend. While satellite remote sensing data (such as GES DISC’s MERRA-2 tavg1_2d_flx_Nx) can provide high-frequency global surface temperature observations, they are still affected by factors such as satellite scanning intervals, cloud cover, instrument failures, and data transmission errors, often resulting in missing data or anomalies. To establish a reliable spatiotemporal temperature field that can support policy-making and scientific research, this study designs a rigorous simulation experiment to compare three time series forecasting methods, SSSD, TSMixer, and RegressionEnsemble, combined with autoFRK for spatial interpolation and prediction, evaluating the differences in performance for data imputation and forecasting.

Research Motivation

The continuous rise in global surface temperatures has become a core issue in climate change research and policy-making. High-resolution and high-temporal-frequency global temperature data are crucial for real-time monitoring of extreme climate events, assessing climate model accuracy, and formulating mitigation and adaptation strategies. However, in practical applications, satellite remote sensing data often suffer from incompleteness or distortion due to observational limitations and technical issues. If left unprocessed, this can affect the reliability of reconstructed temperature fields and future trend forecasts. While existing time series and spatial interpolation methods each have their strengths, there is still a lack of systematic comparative studies on their performance for large-scale spatiotemporal data imputation and forecasting. Therefore, this study aims to conduct a rigorous simulation experiment to assess the applicability and performance differences of various methods in reconstructing global surface temperature data, providing technical references for future climate monitoring and decision-making.

Research Methods

This study is based on the MERRA-2 tavg1_2d_flx_Nx dataset provided by GES DISC, selecting the region with longitude $73^\circ \sim 104^\circ$ and latitude $36^\circ \sim 54^\circ$ as the study area. This area covers northwestern China, western Mongolia, and parts of Kazakhstan, Kyrgyzstan, and Uzbekistan, featuring diverse terrain such as mountains, plateaus, basins, grasslands, lakes, and rivers. The diversity of geographic and climatic conditions results in different yet correlated temperature variation patterns across observation points, serving as the basis for multivariate model inputs. As the study focuses on a relatively small spatiotemporal range, it allows for controlling data heterogeneity while improving model prediction accuracy.

Data Preprocessing

1. Data Extraction and Merging

   • Download the complete 2024 MERRA-2 tavg1_2d_flx_Nx dataset.
   • Merge hourly observation data in chronological order and check the integrity and consistency of time stamps.

2. Missing Value and Outlier Detection

   • Detect missing or anomalous values in the raw data.

3. Data Reshaping

   • Reshape the global dataset into a 3D array of $(24\ \text{hours},\ 366\ \text{days},\ 207{,}936\ \text{locations})$.

4. Experimental Subset Selection

   • Extract the target region (longitude $73^\circ \sim 104^\circ$, latitude $36^\circ \sim 54^\circ$) from the global dataset.
   • Select $(24\ \text{hours},\ 260\ \text{days},\ 1{,}850\ \text{locations})$ as the experimental sample.

5. Defining Known and Unknown Regions

   • Set a fixed pseudo-random seed (seed = 123) to ensure experiment reproducibility.
   • Randomly select 1,500 locations ($\approx$ 81.5%) from the 1,850 locations in the target area as known locations, used for comparing future trend forecasts.
   • Assign the remaining 350 locations ($\approx$ 19.5%) as unknown locations, used for evaluating past spatial imputation performance and future trend forecasts.

The code used is as follows:

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117

# import modules
import os
import numpy as np
import pandas as pd
import xarray as xr
from tqdm import tqdm
import matplotlib.pyplot as plt
import csv
import cartopy.crs as ccrs
import cartopy.feature as cfeature


# functions
def check_data_folder(folder):
    return os.path.exists(folder) and os.path.isdir(folder)

def generate_date_range(start_date, end_date):
    """
    Generate a list of dates from start_date to end_date.
    """
    return pd.date_range(start=start_date, end=end_date, freq='D').strftime('%Y%m%d').tolist()

def load_data(file_path):
    """
    Load data from a NetCDF file.
    """
    if os.path.exists(file_path):
        return xr.open_dataset(file_path)
    else:
        raise FileNotFoundError(f'File not found: {file_path}')


# main program
## check data folder
data_folder = f'..\\..\\..\\surface_air_temperature\\data2024'
if check_data_folder(data_folder):
    print(f'Data folder found: {data_folder}')
else:
    raise FileNotFoundError(f'Data folder not found: {data_folder}')

## load data
start_date = '2024-01-01'
end_date = '2024-12-31'
date_list = generate_date_range(start_date, end_date)

## get location and shape
path = os.path.join(data_folder, f'M2T1NXFLX.5.12.4%3AMERRA2_400.tavg1_2d_flx_Nx.{date_list[0]}.nc4.dap.nc4')
nc4_data = load_data(path)
lat = nc4_data['lat'].values
lon = nc4_data['lon'].values
shape = nc4_data['TLML'].shape
total_locations = shape[1] * shape[2]

## combine data
sample_path = os.path.join(data_folder, f'M2T1NXFLX.5.12.4%3AMERRA2_400.tavg1_2d_flx_Nx.{date_list[0]}.nc4.dap.nc4')
sample_data = load_data(sample_path)
shape_per_file = sample_data['TLML'].shape   # e.g. (24, 361, 576)
time_per_file = len(sample_data['time'])

total_samples = len(date_list)
combined = np.empty((total_samples * shape_per_file[0], *shape_per_file[1:]), dtype=np.float32)
time_list = np.empty(total_samples * time_per_file, dtype=sample_data['time'].dtype)

for i, date in enumerate(tqdm(date_list, desc="Combining")):
    path = os.path.join(data_folder, f'M2T1NXFLX.5.12.4%3AMERRA2_400.tavg1_2d_flx_Nx.{date}.nc4.dap.nc4')
    nc4_data = load_data(path)

    start = i * shape_per_file[0]
    end = (i + 1) * shape_per_file[0]

    combined[start:end] = nc4_data['TLML'].values
    time_list[start:end] = nc4_data['time'].values

print(f'Combined data shape: {combined.shape}')

## reshape data
locations = np.stack(np.meshgrid(lon, lat), axis=-1).reshape(-1, 2)
reshaped_data = combined.reshape(combined.shape[0], -1)
pd.DataFrame(reshaped_data)  # 2d data with time as rows and locations as columns
reshaped_df = pd.DataFrame(reshaped_data, columns=[f"({lon}, {lat})" for lon, lat in locations], index=list(time_list))
time_num = len(time_list)
locations_num = len(locations)

## reshape data to (24, day, location)  (24 hours)
reshaped_df.index = pd.to_datetime(reshaped_df.index)
groups = reshaped_df.groupby(reshaped_df.index.time)
stacked = np.stack([group.to_numpy() for time, group in sorted(groups)])
print(stacked)
time_order = sorted(groups.groups.keys())
print(time_order)

## train set
np.random.seed(123)
valid_mask = (locations[:, 0] >= 73) & (locations[:, 0] <= 104) & \
             (locations[:, 1] >= 36) & (locations[:, 1] <= 54)  # 緯度
valid_indices = np.where(valid_mask)[0]


day_num_train = 250
known_locations_num = valid_indices.shape[0] - 350  # 81.1%
unknown_locations_num = 350                            # 18.9%
locations_num = known_locations_num + unknown_locations_num

locations_index = np.random.choice(valid_indices, size=locations_num, replace=False)
known_locations_index = locations_index[:known_locations_num]
unknown_locations_index = locations_index[known_locations_num:(known_locations_num + unknown_locations_num)]
known_locations_index.sort()
unknown_locations_index.sort()
known_locations_choose = locations[known_locations_index, :]
unknown_locations_choose = locations[unknown_locations_index, :]
stacked_train = stacked[:, :day_num_train, known_locations_index]

future_days = 10
known_real_data = stacked[:, :day_num_train + future_days, known_locations_index]
unknown_real_data = stacked[:, :day_num_train + future_days, unknown_locations_index]

print(f'Load data complete!')

Next, we will generate time series forecasts using the following different models and perform the following operations:

1. Time Series Forecasting

   • Train using the time series data from the known locations with different models.
   • Predict the future 10 days for the known locations, i.e., $24 \times 10 = 240$ time steps.
   • The shape of each prediction is $(24, 10, 1850)$.

2. Spatial Imputation

   • Train the autoFRK model using the time series data and coordinates of the known locations, including the future information predicted by the other models.
   • Impute the full time series values for the unknown locations.

SSSD

SSSD (Structured State Space Diffusion) is a generative time series imputation and forecasting method that combines a conditional diffusion model with a structured state space sequence model (S4). The S4 layers effectively capture long-term dependencies, while the diffusion mechanism enhances the flexibility and expressiveness of imputation, achieving excellent performance in both missing data completion and future trend simulation.

In the diffusion generation mechanism, the imputation task for missing data is formulated as a reverse diffusion process, gradually restoring missing regions from noise. It can handle complex scenarios such as random missing (RM), non-random missing (NRM), and even long blackout missing (BM).

The S4 model is a structured state space model that uses HiPPO theory to initialize the state matrix, enabling efficient and stable capture of long-term structural relationships in time series, thus improving SSSD’s ability to model long-term dependencies and cross-variable correlations.

The SSSD model can handle various types of missing data and forecasting tasks, including:

• Random Missing (RM) Missing points are randomly distributed across multiple time series, with independent missing locations, often used to simulate missing data caused by sporadic observation errors or data transmission failures.

• Random Block Missing (RBM) Missing blocks of random length and position occur in each time series. This can be viewed as a high-dimensional version of non-random missing in a single time series, common in cases such as brief instrument failures or local obstructions.

• Blackout Missing (BM) Similar to random block missing, but all time series have missing data in the same time interval, representing a synchronous missing pattern, often caused by large-scale systematic failures such as satellite outages or data server interruptions.

• Time-Series Forecasting (TF) Similar to blackout missing, but the missing segment is at the end of all time series, corresponding to future periods without observations, used to evaluate the model’s ability to forecast future trends.

Because SSSD requires three-dimensional inputs, here the adjusted dataset shape is $(24, 250, 1850)$. Based on prior experience, using 24-hour data as variables performs significantly better than using 1,850 locations as variables. Therefore, the input shape of the dataset here is adjusted to $(1,850 \text{ locations}, 260 \text{ days}, 24 \text{ hours})$ to achieve more accurate experimental results.

The following are the parameter settings for the SSSD model:

• Model configuration model.yaml In the model settings, both input_channels and output_channels are set to 24, representing 24-dimensional features per input and output (hourly temperatures for 24 hours). The WaveNet part uses 32 residual layers (residual_layers: 32), each with 64 residual channels, and employs skip connections to aggregate outputs from different layers, enhancing multi-scale feature integration.

  For the diffusion model, 200 diffusion steps are set (T: 200), with $\beta$ values linearly increasing from 0.0001 to 0.02, controlling the strength of noise addition and removal.

  The S4 model’s maximum sequence length is set to 250 (s4_max_sequence_length), matching the length of the time series in the dataset, representing 250 days; the state dimension is 64 (s4_state_dim), and both bidirectional computation (s4_bidirectional) and layer normalization (s4_use_layer_norm) are enabled to ensure stability and generalization.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25

wavenet:
  # WaveNet model parameters
  input_channels: 24  # Number of input channels
  output_channels: 24  # Number of output channels
  residual_layers: 32  # Number of residual layers
  residual_channels: 64  # Number of channels in residual blocks
  skip_channels: 64  # Number of channels in skip connections

  # Diffusion step embedding dimensions
  diffusion_step_embed_dim_input: 128  # Input dimension
  diffusion_step_embed_dim_hidden: 512  # Middle dimension
  diffusion_step_embed_dim_output: 512  # Output dimension

  # Structured State Spaces sequence model (S4) configurations
  s4_max_sequence_length: 250  # Maximum sequence length
  s4_state_dim: 64  # State dimension
  s4_dropout: 0.0  # Dropout rate
  s4_bidirectional: true  # Whether to use bidirectional layers
  s4_use_layer_norm: true  # Whether to use layer normalization

diffusion:
  # Diffusion model parameters
  T: 200  # Number of diffusion steps
  beta_0: 0.0001  # Initial beta value
  beta_T: 0.02  # Final beta value

• Model Training Configuration training.yaml

  In the training settings, the batch size (batch_size) is set to 300. Since there are 1,500 known locations in the training data, each iteration is divided into $1500 \div 300 = 5$ batches to ensure that each iteration covers the entire training set. The maximum number of iterations (n_iters) is set to 3,800, and the learning rate (learning_rate) is 0.0005.

  The training strategy uses only generating missing values (only_generate_missing: true), and applies a time series forecasting missing mask (masking: "forecast"), which treats the end of the sequence as unknown for prediction. The length of each masked segment (missing_k) is 10.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18

# Training configuration
batch_size: 300  # Batch size
output_directory: "./results/surface air temperature/control"  # Output directory for checkpoints and logs
ckpt_iter: "max"  # Checkpoint mode (max or min)
iters_per_ckpt: 100  # Checkpoint frequency (number of epochs)
iters_per_logging: 100  # Log frequency (number of iterations)
n_iters: 3800  # Maximum number of iterations
learning_rate: 0.0005  # Learning rate

# Additional training settings
only_generate_missing: true  # Generate missing values only
use_model: 2  # Model to use for training
masking: "forecast"  # Masking strategy for missing values
missing_k: 10  # Number of missing values

# Data paths
data:
  train_path: "./datasets/surface air temperature/train/train.npy"  # Path to training data

• Model Imputation Configuration inference.yaml

  The inference settings are similar to the training configuration, with a batch size of 300 (batch_size), and the model is loaded from the best checkpoint path (ckpt_path) saved during training. The number of inference repetitions (trials) is set to 1, and each time generation is performed only on the missing parts. The missing mask strategy and length remain consistent with training to ensure uniformity.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15

# Inference configuration
batch_size: 300  # Batch size for inference
output_directory: "./results/surface air temperature/inference/control"  # Output directory for inference results
ckpt_path: "./results/surface air temperature/control"  # Path to checkpoint for inference
trials: 1 # Replications

# Additional training settings
only_generate_missing: true  # Generate missing values only
use_model: 2  # Model to use for training
masking: "forecast"  # Masking strategy for missing values
missing_k: 10  # Number of missing values

# Data paths
data:
  test_path: "./datasets/surface air temperature/test/test.npy"  # Path to test data

TSMixer

Unlike deep generative models like SSSD, TSMixer (Time-Series Mixer) is a time series forecasting model based on a multilayer perceptron (MLP) architecture, designed to efficiently capture temporal and feature dimension correlations within time series for multivariate forecasting. TSMixer models time series data through stacked MLP mixer layers that alternate fusion along the time and feature dimensions.

The core of TSMixer consists of multiple mixer layers, each first mixing along the time dimension and then along the feature dimension. This enables the model to effectively capture temporal dependencies and feature correlations within the time series. The entire model is composed solely of MLPs, avoiding computational bottlenecks present in RNNs or attention mechanisms, thereby improving training efficiency and prediction speed.

Additionally, TSMixer supports various types of auxiliary variables, including past covariates, future covariates, and static covariates, providing high flexibility for multivariate and complex forecasting tasks. In several long-term forecasting benchmarks, TSMixer outperforms traditional Transformer-based models in both computational efficiency and accuracy.

TSMixer is suitable for:

• Multivariate time series forecasting Able to handle correlations among multiple variables, applicable to fields such as weather forecasting and financial market analysis.

• Long-term forecasting tasks Provides stable and accurate predictions in long-term forecasts.

• Resource-constrained environments Due to its efficient computational performance, TSMixer is well-suited for scenarios with limited computing resources.

Below, we use the PyTorch-based darts module’s TSMixerModel to train and forecast multivariate time series on a daily (24-hour) basis.

In this code, input_chunk_length=30 means the model observes the past 30 days of data, while output_chunk_length=10 predicts the next 10 days. The training runs for 3,800 iterations.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35

from darts.models import TSMixerModel
from darts import TimeSeries
from datetime import datetime
from tqdm import tqdm
real_data = known_real_data
location_choose = known_locations_choose
locations_index = known_locations_index
inference = np.array([[[0.0] * real_data.shape[2]] * future_days] * real_data.shape[0])
start_time = datetime.now()

for i in tqdm(range(real_data.shape[0])):
    train_set = real_data[i, :day_num_train, :]
    train_set = pd.DataFrame(train_set, index=pd.to_datetime(date_list[:day_num_train]), columns=[f"loc_{i}" for i in range(known_locations_num)])
    train_set = TimeSeries.from_dataframe(train_set)

    test_set = real_data[i, day_num_train:, :]
    test_set = pd.DataFrame(test_set, index=pd.to_datetime(date_list[day_num_train:day_num_train + future_days]), columns=[f"loc_{i}" for i in range(known_locations_num)])
    test_set = TimeSeries.from_dataframe(test_set)

    model = TSMixerModel(
        input_chunk_length=30,
        output_chunk_length=10,
        n_epochs=3800,
        dropout=0.0005,
        use_reversible_instance_norm=True,
        random_state=42,
        pl_trainer_kwargs={"accelerator": "gpu"}
    )

    model.fit(train_set)

    forecast = model.predict(future_days)
    inference[i] = forecast.values()

print(f'Inference complete! Time taken: {datetime.now() - start_time}')

RegressionEnsemble

RegressionEnsemble is an ensemble forecasting model in the darts module that uses regression models (such as linear regression) to combine outputs from multiple base forecasting models, improving overall prediction accuracy. This model employs a stacking technique, using the predictions from multiple base models as features, and trains a regression model to learn the optimal fusion weights.

The model can simultaneously forecast multiple variables and integrate various auxiliary information, including past known variables, future known variables, and static covariates, enhancing prediction accuracy. RegressionEnsemble is suitable for:

• Multivariate time series forecasting Effectively combines predictions from different models when forecasting multiple variables simultaneously.

• Long-term forecasting tasks Improves prediction stability and accuracy by aggregating multiple model results.

• Resource-constrained environments Due to the computational efficiency of regression fusion, it is particularly suitable for environments with limited computing resources.

Below, we use the PyTorch-based darts module’s RegressionEnsembleModel to train and forecast multivariate time series on a daily (24-hour) basis.

In this code, three base forecasting models are used to build the regression ensemble model. First, NaiveSeasonal(K=7) is a simple seasonal model based on the assumption that the time series exhibits similar patterns every 7 days. Second, LinearRegressionModel(lags=30) is a linear regression model using data from the past 30 time points as features to capture trends and changes in the time series. Lastly, the NaiveDrift() model predicts simple trends based on the overall drift of the sequence.

The predictions from these three base models are used as input features for the regression ensemble to learn how to optimally combine them, enhancing overall forecast accuracy. The parameter regression_train_n_points=30 specifies that the fusion model uses predictions from the past 30 time points to train the regression model, enabling it to more effectively capture temporal patterns and generate accurate forecasts.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36

from darts.models import RegressionEnsembleModel, NaiveSeasonal, LinearRegressionModel, NaiveDrift
from darts import TimeSeries
from datetime import datetime
from tqdm import tqdm
real_data = known_real_data
location_choose = known_locations_choose
locations_index = known_locations_index
inference = np.array([[[0.0] * real_data.shape[2]] * future_days] * real_data.shape[0])
start_time = datetime.now()

for i in tqdm(range(real_data.shape[0])):
    train_set = real_data[i, :day_num_train, :]
    train_set = pd.DataFrame(train_set, index=pd.to_datetime(date_list[:day_num_train]), columns=[f"loc_{i}" for i in range(known_locations_num)])
    train_set = TimeSeries.from_dataframe(train_set)

    test_set = real_data[i, day_num_train:, :]
    test_set = pd.DataFrame(test_set, index=pd.to_datetime(date_list[day_num_train:day_num_train + future_days]), columns=[f"loc_{i}" for i in range(known_locations_num)])
    test_set = TimeSeries.from_dataframe(test_set)

    base_models = [
        NaiveSeasonal(K=7),
        LinearRegressionModel(lags=30),
        NaiveDrift()
    ]

    model = RegressionEnsembleModel(
        forecasting_models=base_models,
        regression_train_n_points=30
    )

    model.fit(train_set)

    forecast = model.predict(n=future_days)
    inference[i] = forecast.values()

print(f'Inference complete! Time taken: {datetime.now() - start_time}')

autoFRK

autoFRK (Automatic Fixed Rank Kriging) is an efficient spatial interpolation and forecasting method based on spatial statistical theory. It combines the dimension reduction technique of fixed rank kriging (FRK), using multi-scale basis functions to capture different spatial variation features in data, and automatically selects model parameters.

Mathematically, the autoFRK model can be expressed as:

$$ z[t] = \mu + G \cdot w[t] + \eta[t] + e[t], \quad w[t] \sim N(0, M), \quad e[t] \sim N(0, s \cdot D); \quad t = 1, \cdots, T, $$

where $z[t]$ is the observed (partial) data vector at $n$ locations; $\mu$ is a constant mean vector of length $n$; $D$ is a known $n \times n$ matrix; $G$ is a known $n \times K$ matrix; $\eta[t]$ is a random vector of length $n$ corresponding to a spatially stationary process; and $w[t]$ is an unobserved random weight vector of length $K$.

Parameters are estimated via maximum likelihood with closed-form expressions. The basis function matrix $G$ is constructed using ordered thin-plate spline functions, and the number of bases is selected by Akaike’s information criterion (AIC).

autoFRK represents the spatial random field via basis functions $G$, compressing high-dimensional spatial data into a low-dimensional coefficient space, greatly reducing computational complexity. It is especially suitable for large-scale spatial data analysis. By maximizing the likelihood function, autoFRK effectively captures spatial structure and provides accurate spatial predictions along with uncertainty quantification.

In this study, autoFRK is used to perform spatial interpolation and imputation on temporal prediction results generated by time series models (such as SSSD, TSMixer, RegressionEnsemble). This spatiotemporal integration helps fill data gaps and enhances the completeness and accuracy of global surface temperature monitoring.

The application procedure of autoFRK in this study is as follows:

1. Input the point-wise temporal predictions generated by time series models.
2. Perform spatial basis function expansion and parameter estimation via autoFRK.
3. Conduct spatial interpolation and prediction for unknown locations to fill missing data areas.
4. Combine temporal and spatial prediction results to provide a complete global surface temperature spatiotemporal field.

The code for spatial imputation is as follows:

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129

# load data
library(reticulate)
np <- import("numpy")

data_path = "../train_frk.npy"
data = np$load(data_path)

real_path = "../real_data.npy"
real_data = np$load(real_path)
known_data = real_data[ , , 1:1500]
unknown_data = real_data[ , , 1501:1850]

known_loc_path = "../known_location.npy"
known_locs = np$load(known_loc_path)

unknown_loc_path = "../unknown_location.npy"
unknown_locs = np$load(unknown_loc_path)

# autoFRK
library(autoFRK)
library(dplyr)

start_time <- Sys.time()

unknown_inference_shape = c(dim(data)[1:2],  dim(known_locs)[1] + dim(unknown_locs)[1])
unknown_inference <- array(NA, dim = unknown_inference_shape)
mrts_basis <- NULL

locs = rbind(known_locs, unknown_locs)
total_iter <- 24 * dim(data)[2]
pb <- txtProgressBar(min = 0, max = total_iter, style = 3)
iter <- 0

for (time_ in 1:24) {
  hour_data <- data[time_, , ]
  
  for (day_ in 1:dim(hour_data)[1]) {
    iter <- iter + 1
    iter_start <- Sys.time()
    
    day_data = hour_data[day_, ]
    
    if (is.null(mrts_basis)) {
      model = autoFRK(data = day_data, loc = known_locs)
      mrts_basis = model$G
    } else {
      model = autoFRK(data = day_data, loc = known_locs, G = mrts_basis)
    }
    
    pred = predict.FRK(object = model, newloc = locs)
    unknown_inference[time_, day_, ] = pred$pred.value
    
    elapsed <- as.numeric(difftime(Sys.time(), start_time, units = "secs"))
    avg_time <- elapsed / iter
    remaining <- avg_time * (total_iter - iter)
    est_done <- format(Sys.time() + remaining, "%H:%M:%S")
    
    setTxtProgressBar(pb, iter)
    cat(sprintf(" | Est. done at: %s | Remaining: %ds\r", est_done, round(remaining)))
  }
}

close(pb)

end_time <- Sys.time()
cat("\n")
cat("Total time elapsed:", format(end_time - start_time))
cat("\n")

# result
y_inf <- unknown_inference
real <- real_data
train <- known_data
test <- unknown_data

mspe <- function(pred, true) mean((pred - true)^2)
rmspe <- function(pred, true) sqrt(mean((pred - true)^2))
mape <- function(pred, true) mean(abs(pred - true))

mspe_p <- function(pred, true) mean((pred - true)^2 / true)
rmspe_p <- function(pred, true) sqrt(mean((pred - true)^2 / true))
mape_p <- function(pred, true) mean(abs(pred - true) / true)

future_days = 10
future_idx <- (dim(y_inf)[2] - future_days + 1):dim(y_inf)[2]
past_idx <- 1:(dim(y_inf)[2] - future_days)
known_idx <- 1:dim(known_locs)[1]
unknown_idx <- (dim(known_locs)[1] + 1):(dim(known_locs)[1] + dim(unknown_locs)[1])

compute_metrics <- function(pred, true) {
  c(
    MSPE = mspe(pred, true),
    RMSPE = rmspe(pred, true),
    `MSPE%` = mspe_p(pred, true),
    `RMSPE%` = rmspe_p(pred, true),
    MAPE = mape(pred, true),
    `MAPE%` = mape_p(pred, true)
  )
}

result_table <- data.frame(
  row.names = c('MSPE', 'RMSPE', 'MSPE%', 'RMSPE%', 'MAPE', 'MAPE%'),

  `ALL Locs & All Time` = compute_metrics(y_inf, real),
  `Known Locs & All Time` = compute_metrics(y_inf[, , known_idx], train),
  `Unknown Locs & All Time` = compute_metrics(y_inf[, , unknown_idx], test),

  `ALL Locs & Future` = compute_metrics(y_inf[, future_idx, ], real[, future_idx, ]),
  `Known Locs & Future` = compute_metrics(y_inf[, future_idx, known_idx], train[, future_idx, ]),
  `Unknown Locs & Future` = compute_metrics(y_inf[, future_idx, unknown_idx], test[, future_idx, ]),

  `ALL Locs & Past` = compute_metrics(y_inf[, past_idx, ], real[, past_idx, ]),
  `Known Locs & Past` = compute_metrics(y_inf[, past_idx, known_idx], train[, past_idx, ]),
  `Unknown Locs & Past` = compute_metrics(y_inf[, past_idx, unknown_idx], test[, past_idx, ])
)

print(result_table)


# save to .npy

# combined data
output_matrix = unknown_inference
output_matrix %>% dim() %>% print()

# save to .npy
output_matrix = output_matrix %>% r_to_py()
save_path = "../plot_frk.npy"
np$save(save_path, output_matrix)

Experimental Results

In the experiments above, except for the RegressionEnsemble model, all models were trained for 3,800 iterations. The data format for all three models was the same; however, since TSMixer and RegressionEnsemble only accept two-dimensional data as input, training and forecasting were performed using a loop.

SSSD + autoFRK

For conclusions, please refer to the section 20250731 meeting Experiment 3.

TSMixer + autoFRK

For conclusions, please refer to the section 20250808 meeting TSMixerModel + autoFRK.

RegressionEnsemble + autoFRK

For conclusions, please refer to the section 20250808 meeting RegressionEnsembleModel + autoFRK.

Conclusion

From the above experiments, it can be observed that SSSD + autoFRK demonstrates relatively stable performance in time series forecasting (Future), especially in the unknown locations’ future predictions (Unknown Locs & Future), where metrics such as MSPE, RMSPE, and MAPE are significantly lower than those of TSMixer + autoFRK, and the gap compared to known locations (Known Locs) is minimal. This indicates better generalization ability at unobserved locations. However, the SSSD model incurs higher time costs, making it less suitable for scenarios requiring real-time forecasting. Increasing the number of iterations without significantly raising computational burden may further enhance its predictive power.

In contrast, TSMixer + autoFRK shows notably higher MSPE and RMSPE in time series forecasting (e.g., MSPE as high as 75.99), especially larger errors in unknown regions, indicating substantial prediction bias at unobserved locations. This model may require adjustments to improve generalization ability, such as increasing iteration count, refining training strategies, or incorporating additional features to close the gap with other methods.

The overall performance of RegressionEnsemble + autoFRK lies between SSSD and TSMixer. Although MSPE and RMSPE in unknown areas are higher than SSSD, they are much lower than TSMixer, indicating moderate to good generalization ability. Moreover, RegressionEnsemble has significantly shorter inference time than SSSD. Introducing more efficient base forecasting models without significantly increasing computational costs could potentially allow it to outperform SSSD in predictive accuracy.

Epilogue

Environment

• Local Operating System: Windows 11 24H2
  • Programming Language: Python 3.12.9
• Computing Platform: National Center for High-Performance Computing (NCHC) – Taiwan AI Cloud
  • Operating System: Ubuntu
  • Miniconda
  • GPU: NVIDIA Tesla V100 32GB GPU
  • CUDA 12.8 driver
  • Programming Language: Python 3.10.16 for Linux

Further Learning

• I tested this project on the GitHub repository.

References

• Global Modeling and Assimilation Office (GMAO). (2015). MERRA-2 tavg1_2d_flx_Nx: 2d,1-Hourly,Time-Averaged,Single-Level,Assimilation,Surface Flux Diagnostics (Version 5.12.4) [Dataset]. Goddard Earth Sciences Data and Information Services Center (GES DISC). Retrieved from https://doi.org/10.5067/7MCPBJ41Y0K6

• Juan Lopez Alcaraz, Nils Strodthoff (2022). Diffusion-based time series imputation and forecasting with structured state space models. Transactions on Machine Learning Research. Retrieved from https://openreview.net/forum?id=hHiIbk7ApW

• SSSD (2022). GitHub. Retrieved from https://github.com/AI4HealthUOL/SSSD

• SSSD_CP (2024). GitHub. Retrieved from https://github.com/egpivo/SSSD_CP

• Unit8 SA (n.d.). Time Series Made Easy in Python. Darts. Retrieved from https://unit8co.github.io/darts/index.html

• darts (2025). GitHub. Retrieved from https://github.com/unit8co/darts

• Tzeng, S., & Huang, H. C. (2018). Resolution Adaptive Fixed Rank Kriging. Technometrics, 60(2), 198–208. Retrieved from https://doi.org/10.1080/00401706.2017.1345701

• autoFRK (2024). GitHub. Retrieved from https://github.com/egpivo/autoFRK

• Si-An Chen, Chun-Liang Li, Nate Yoder, Sercan O. Arik, and Tomas Pfister. (2023). TSMixer: An all-MLP architecture for time series forecasting. arXiv. Retrieved from https://arxiv.org/abs/2303.06053