Markov Chain Gaussian Fields Simulation and Parameter Estimation.
mcgf
The goal of mcgf is to provide easy-to-use functions for simulating and fitting covariance models. It provides functions for simulating (regime-switching) Markov chain Gaussian fields with covariance functions of the Gneiting class by simple kriging. Parameter estimation methods such as weighted least squares and maximum likelihood estimation are available. Below is an example of simulating and estimation parameters for an MCGF.
Installation
You can install the development version of mcgf from GitHub with:
# install.packages("devtools")
devtools::install_github("tianxia-jia/mcgf")
Data Simulation
To simulate an MCGF with fully symmetric covariance structure, we begin with simulating 10 locations randomly.
library(mcgf)
set.seed(123)
h <- rdists(10)
Next, we simulate an MCGF with the general stationary covariance structure. In this example the covariance structure is a convex combination of a base separable model and a Lagrangian model account for asymmetry.
N <- 1000
lag <- 5
par_base <- list(
par_s = list(nugget = 0, c = 0.001, gamma = 0.5),
par_t = list(a = 0.5, alpha = 0.8)
)
par_lagr <- list(v1 = 200, v2 = 200, k = 2)
sim1 <- mcgf_sim(
N = N,
base = "sep",
lagrangian = "lagr_tri",
par_base = par_base,
par_lagr = par_lagr,
lambda = 0.2,
dists = h,
lag = lag
)
sim1 <- sim1[-c(1:(lag + 1)), ]
rownames(sim1) <- 1:nrow(sim1)
sim1 <- list(data = sim1, dists = h)
Parameter Estimation
Create an mcgf object
To estimate parameters, we need to calculate auto-correlations and cross-correlations. Let’s first create an mcgf object. The mcgf class extends the data.frame with more attributes.
sim1_mcgf <- mcgf(sim1$data, dists = sim1$dists)
#> `time` is not provided, assuming rows are equally spaced temporally.
Then the acfs and ccfs can be added to this object as follows.
sim1_mcgf <- add_acfs(sim1_mcgf, lag_max = lag)
sim1_mcgf <- add_ccfs(sim1_mcgf, lag_max = lag)
Estimate base model
To perform parameter estimation, we can start with estimating the parameters for spatial and temporal models.
fit_spatial <- fit_base(
sim1_mcgf,
model = "spatial",
lag = lag,
par_init = c(c = 0.001, gamma = 0.5),
par_fixed = c(nugget = 0)
)
fit_spatial$fit
#> $par
#> c gamma
#> 0.001160802 0.500000000
#>
#> $objective
#> [1] 1.640593
#>
#> $convergence
#> [1] 0
#>
#> $iterations
#> [1] 8
#>
#> $evaluations
#> function gradient
#> 21 20
#>
#> $message
#> [1] "both X-convergence and relative convergence (5)"
fit_temporal <- fit_base(
sim1_mcgf,
model = "temporal",
lag = lag,
par_init = c(a = 0.3, alpha = 0.5)
)
fit_temporal$fit
#> $par
#> a alpha
#> 0.6528906 0.7560970
#>
#> $objective
#> [1] 0.004306706
#>
#> $convergence
#> [1] 0
#>
#> $iterations
#> [1] 18
#>
#> $evaluations
#> function gradient
#> 23 43
#>
#> $message
#> [1] "relative convergence (4)"
Alternatively, we can fit the separable model all at once:
fit_sep <- fit_base(
sim1_mcgf,
model = "sep",
lag = lag,
par_init = c(
c = 0.001,
gamma = 0.5,
a = 0.5,
alpha = 0.5
),
par_fixed = c(nugget = 0)
)
fit_sep$fit
#> $par
#> c gamma a alpha
#> 0.001154864 0.500000000 0.624551338 0.735490605
#>
#> $objective
#> [1] 3.488305
#>
#> $convergence
#> [1] 0
#>
#> $iterations
#> [1] 18
#>
#> $evaluations
#> function gradient
#> 49 88
#>
#> $message
#> [1] "relative convergence (4)"
we can also estimate the parameters using MLE:
fit_sep2 <- fit_base(
sim1_mcgf,
model = "sep",
lag = lag,
par_init = c(
c = 0.001,
gamma = 0.5,
a = 0.5,
alpha = 0.5
),
par_fixed = c(nugget = 0),
method = "mle",
)
fit_sep2$fit
#> $par
#> c gamma a alpha
#> 0.001197799 0.500000000 0.804621207 1.000000000
#>
#> $objective
#> [1] -11520.04
#>
#> $convergence
#> [1] 0
#>
#> $iterations
#> [1] 17
#>
#> $evaluations
#> function gradient
#> 55 78
#>
#> $message
#> [1] "relative convergence (4)"
Now we will add the base model to x_mcgf:
sim1_mcgf <- add_base(sim1_mcgf, fit_base = fit_sep)
To print the current model, we do
model(sim1_mcgf)
#> ----------------------------------------
#> Model
#> ----------------------------------------
#> - Time lag: 5
#> - Scale of time lag: 1
#> - Forecast horizon: 1
#> ----------------------------------------
#> Base
#> ----------------------------------------
#> - Base model: sep
#> - Parameters:
#> c gamma a alpha nugget
#> 0.001154864 0.500000000 0.624551338 0.735490605 0.000000000
#>
#> - Fixed parameters:
#> nugget
#> 0
#>
#> - Parameter estimation method: wls
#>
#> - Optimization function: nlminb
#> ----------------------------------------
#> Lagrangian
#> ----------------------------------------
#> - Lagrangian model:
#> - Parameters:
#> NULL
#>
#> - Fixed parameters:
#> NULL
#>
#> - Parameter estimation method:
#>
#> - Optimization function:
Estimate the Lagrangian model
Similarly, we can estimate the parameters for the Lagrangian component by
fit_lagr <- fit_lagr(
sim1_mcgf,
model = "lagr_tri",
par_init = c(v1 = 300, v2 = 300, lambda = 0.15),
par_fixed = c(k = 2)
)
fit_lagr$fit
#> $par
#> lambda v1 v2
#> 0.1757035 232.0852117 203.8869305
#>
#> $objective
#> [1] 1.627017
#>
#> $convergence
#> [1] 0
#>
#> $iterations
#> [1] 32
#>
#> $evaluations
#> function gradient
#> 35 126
#>
#> $message
#> [1] "relative convergence (4)"
We can add the Lagrangian model by
sim1_mcgf <- add_lagr(sim1_mcgf, fit_lagr = fit_lagr)
Finally we may print the final model:
model(sim1_mcgf)
#> ----------------------------------------
#> Model
#> ----------------------------------------
#> - Time lag: 5
#> - Scale of time lag: 1
#> - Forecast horizon: 1
#> ----------------------------------------
#> Base
#> ----------------------------------------
#> - Base model: sep
#> - Parameters:
#> c gamma a alpha nugget
#> 0.001154864 0.500000000 0.624551338 0.735490605 0.000000000
#>
#> - Fixed parameters:
#> nugget
#> 0
#>
#> - Parameter estimation method: wls
#>
#> - Optimization function: nlminb
#> ----------------------------------------
#> Lagrangian
#> ----------------------------------------
#> - Lagrangian model: lagr_tri
#> - Parameters:
#> lambda v1 v2 k
#> 0.1757035 232.0852117 203.8869305 2.0000000
#>
#> - Fixed parameters:
#> k
#> 2
#>
#> - Parameter estimation method: wls
#>
#> - Optimization function: nlminb
Kriging forecast
This package provides kriging forecasts (and intervals) for empirical, base, and general stationary models.
# Empirical model
fit_emp <-
krige(sim1_mcgf,
model = "empirical",
interval = TRUE
)
rmse_emp <- sqrt(mean(colMeans((sim1_mcgf - fit_emp$fit)^2, na.rm = T)))
# Base separable model
fit_base <-
krige(sim1_mcgf,
model = "base",
interval = TRUE
)
rmse_base <-
sqrt(mean(colMeans((sim1_mcgf - fit_base$fit)^2, na.rm = T)))
# Stationary model
fit_stat <-
krige(sim1_mcgf,
model = "all",
interval = TRUE
)
rmse_stat <-
sqrt(mean(colMeans((sim1_mcgf - fit_stat$fit)^2, na.rm = T)))
rmse <- c(rmse_emp, rmse_base, rmse_stat)
names(rmse) <- c("Empirical", "Separable", "Stationary")
rmse
#> Empirical Separable Stationary
#> 0.7212175 0.7685016 0.7355458