Spatiotemporal stochastic volatility model -- Environmental risk modelling

In this live script, we demonstrate the usage of the dynamic stochastic volatility model using an empirical example from environmental science. As for applications in financial economics, the volatility of the process can be interpreted as risks (i.e., environmental risks in our case).

The full paper can be found on arXiv: arXiv:2211.03178

Load packages, functions and Agrimonia dataset

We analyze the variation of fine dust concentrations of particles having a diameter less than 10 micrometers, PM10, in Lombardy, Northern Italy. Bounded by the Alps to the west, the region experiences reduced wind circulation, contributing to Lombardy being among the European regions with the lowest air quality. We use the daily PM10 concentrations from 1.1.2021 to 31.12.2021 from the official monitoring stations of the regional environmental authority, ARPA Lombardia. The data are open-source provided by the Agrimonia project: https://zenodo.org/record/7956006 (10.5281/zenodo.7956006).

clear all;
 
addpath(genpath("toolbox_may2021"))
 
load("AgrImOnIA_v3.0.0/AgrImOnIA_Dataset_v_3_0_0.mat")
agrimonia = AgrImOnIADataset;

Reshaping of data, imputation and standardization

time = agrimonia.Time;
t = length(unique(time));
n = length(unique(agrimonia.IDStations));
locations = unique(agrimonia.IDStations);
 
coords = zeros(n, 3);
 
for i = 1:length(locations)
    idx = find(agrimonia.IDStations == locations(i));
    coords(i, 1) = i;
    coords(i, 2) = agrimonia.Latitude(idx(1));
    coords(i, 3) = agrimonia.Longitude(idx(1));
end
 
ymat = reshape(agrimonia.AQ_pm10, t, n)';
 
t2 = 2192;
t1 = 1827;
ymat = ymat(:, (t1+1):t2);
n2ind = find(mean(isnan(ymat), 2) ~= 1);
n2 = size(n2ind, 1);
ymat = ymat(n2ind, :);
 
ymat(isnan(ymat)) = nanmean(ymat, "all"); % na imputation
ymat(ymat == 0) = min(ymat(ymat > 0));
 
y = reshape(ymat, n2*(t2-t1), []); % reshape
 
covmat_blhmax = (reshape(agrimonia.WE_blh_layer_max, t, n)' - mean(agrimonia.WE_blh_layer_max, "all")) / sqrt(var(agrimonia.WE_blh_layer_max));
covmat_temp   = (reshape(agrimonia.WE_temp_2m, t, n)' - mean(agrimonia.WE_temp_2m, "all")) / sqrt(var(agrimonia.WE_temp_2m));
covmat_humid  = (reshape(agrimonia.WE_rh_mean, t, n)' - mean(agrimonia.WE_rh_mean, "all")) / sqrt(var(agrimonia.WE_rh_mean));
covmat_press  = (reshape(agrimonia.WE_surface_pressure, t, n)' - mean(agrimonia.WE_surface_pressure, "all")) / sqrt(var(agrimonia.WE_surface_pressure));
 
x = [ reshape(covmat_blhmax, n*t, []) reshape(covmat_temp, n*t, []) reshape(covmat_humid, n*t, []) reshape(covmat_press, n*t, [])]; % nT x p-dimensional matrix of regressors
xmat = reshape(x, n, t, 4);
xmat = xmat(n2ind, (t1+1):t2, :);
x = reshape(xmat, n2*(t2-t1), 4);
 
distmatfull = distance_miles(coords(:, 2), coords(:, 3));
coordsfull = coords;
 
coords = coordsfull(n2ind, :);
distmat = distance_miles(coords(:, 2), coords(:, 3));

Descriptive plots

figure;

time_idx = string(unique(time));

ts1 = timeseries(median(ymat, 1), char(time_idx((t1+1):t2)));

ts2 = timeseries(quantile(ymat, 0.05, 1), char(time_idx((t1+1):t2)));

ts3 = timeseries(quantile(ymat, 0.95, 1), char(time_idx((t1+1):t2)));

ts1.Name = 'PM10 concentrations';

ts1.TimeInfo.Units = "days";

ts1.TimeInfo.StartDate = '2021-01-01';

ts1.TimeInfo.Format = 'yyyy-mm-dd';

ts2.TimeInfo.Units = "days";

ts2.TimeInfo.StartDate = '2021-01-01';

ts2.TimeInfo.Format = 'yyyy-mm-dd';

ts3.TimeInfo.Units = "days";

ts3.TimeInfo.StartDate = '2021-01-01';

ts3.TimeInfo.Format = 'yyyy-mm-dd';

p1 = plot(ts1, "-", 'LineWidth', 2, 'Color', [0,0,0.6]); hold on;

plot(ts2, "--", 'LineWidth', 1, 'Color', [0,0,0.6]); hold on;

plot(ts3, "--", 'LineWidth', 1, 'Color', [0,0,0.6]); hold on;

grid on;

hold off;

title('Median PM10 Concentrations','FontSize', 14, 'interpreter', 'latex');

figure;

p2 = geoscatter(coords(:,2)', coords(:,3)', 50, median(ymat, 2), "filled");

geobasemap 'street';

colormap(turbo);

colorbar

hold off;

title('Median PM10 Concentrations','FontSize', 14, 'interpreter', 'latex');