Obtaining covariance matrix from space time variogram

Hi all,

I'm trying to obtain a full covariance matrix for a STIDF object
(spacetime library).
I use the vgmST function from gstat to specify a functional form based
on regression residuals of a subset.

While the krigeST function provides the option to compute such a full
covariance matrix, this only works for global kriging.
As this requires a lot of computing power, this is no option for me.
Hence, I am looking for a method to obtain the covariance matrix
directly.

While I found various functions to do so with purely spatial
variograms, this was not the case for space time variograms.
Any help will be appreciated.

Kind regards,
Philibert

Example Script:

#--R-Script--#

library(sp)
library(spacetime)
library(gstat)
data(air)

suppressWarnings(proj4string(stations) <- CRS(proj4string(stations)))
rural = STFDF(stations, dates, data.frame(PM10 = as.vector(air)))
rr <- rural[,"2005-06-01/2005-09-01"]
rr <- as(rr,"STSDF")

vgram <- variogramST(formula = PM10 ~ 1, data = rr, tunit = "days",
tlags = seq(0,3,1), cutoff = 500, cores = 7)

sumMetricVgm <- vgmST("sumMetric",
                       space = vgm( 8.5, "Exp", 196.6,  3),
                       time  = vgm( 25, "Sph",   1.1,  2),
                       joint = vgm(5, "Exp", 136.6, 12),
                       stAni = 51.7)

VAR <- fit.StVariogram(vgram, sumMetricVgm,
                        lower = c(0,0,0,0,0,0,0,0,0,0),,
                        method="L-BFGS-B",
                        fit.method = 7)
plot(vgram,VAR, all = TRUE)

#Thus far I have tried the krige method, but this is not an option as
the newdata contains a larger set.
kr <- krigeST(PM10 ~ 1, data = rr, newdata = rr, VAR, computeVar =
TRUE, fullCovariance = TRUE)

#While the local kriging option is computably feasable it does not
provide full covariance matrices.
lkr <- krigeST(PM10 ~ 1, data = rr, newdata = rural, VAR, nmax = 10,
computeVar = TRUE, fullCovariance = TRUE)

#--R-Script--#