This vignette describes the use of the DCSmooth-package and its functions. The DCSmooth package provides some tools for non- and semiparametric estimation of the mean surface m of an observed sample of some function y(x, t) = m(x, t) + ε(x, t).
This package is accompanying Schäfer and Feng (2021). Boundary modification procedures for local polynomial smoothing are considered by Feng and Schäfer (2021). For a more comprehensive and detailed description of the algorithms than in section 4, please refer to these papers.
The DCSmooth contains the following functions, methods and data sets:
| Functions | |
|---|---|
set.options() |
Define options for the
dcs()-function. |
dcs() |
Nonparametric estimation of the expectation function of
a matrix Y. Includes automatic iterative plug-in bandwidth
selection. |
surface.dcs() |
3d-plot for the surfaces of an
"dcs"-object. |
sarma.sim() |
Simulate a SARMA-model. |
sarma.est() |
Estimate the parameters of a SARMA-model. |
sfarima.sim() |
Simulate a SFARIMA-model. |
sfarima.est() |
Estimate the parameters of a SFARIMA-model. |
kernel.assign() |
Assign a pointer to a kernel function. |
kernel.list() |
Print list of kernels available in the package. |
| Methods/Generics | |
|---|---|
summary.dcs() |
Summary statistics for an object of class
"dcs". |
print.dcs() |
Print an object of class "dcs". |
plot.dcs() |
Plot method for an "dcs"-object, returns
contour plot. |
residuals.dcs() |
Returns the residuals of the regression from an
"dcs"-object. |
print.summary_dcs() |
Print an object of class "summary_dcs",
which inherits from summary.dcs(). |
print.set_options() |
Prints an object of class "dcs_options",
which inherits from set.options() |
summary.sarma() |
Summary statistics for an object of class
"sarma" |
print.summary_sarma() |
Prints an object of class "summary_sarma",
which inherits from summary.sarma() |
summary.sfarima() |
Summary statistics for an object of class
"sfarima" |
print.summary_sfarima() |
Prints an object of class
"summary_sfarima", which inherits from
summary.sfarima() |
| Data | |
|---|---|
y.norm1 |
A surface with a single gaussian peak. |
y.norm2 |
A surface with two gaussian peaks. |
y.norm3 |
A surface with two gaussian ridges. |
temp.nunn |
Temperatures in Nunn, CO observed in 2020 in 5 minute intervals. (Source: NOAA) |
temp.yuma |
Temperatures in Yuma, AZ observed in 2020 in 5 minute intervals. (Source: NOAA) |
wind.nunn |
Wind speed in Nunn, CO observed in 2020 in 5 minute intervals. (Source: NOAA) |
wind.yuma |
Wind speed in Yuma, AZ observed in 2020 in 5 minute intervals. (Source: NOAA) |
returns.alv |
5 minute returns of Allianz SE from 2007 to 2010 |
volumes.alv |
5 minute volumes of Allianz SE from 2007 to 2010 |
set.options()This auxiliary function is used to set the options for the
dcs function. An object of class dcs_options
is created and should be used as dcs_options- argument in
the dcs function.
Arguments of set.options() are
"KR") and local polynomial
regression ("LP"), which is the default value.M, MW
or T. The value k
is the kernel order, μ is the
smoothness degree and ν the
derivative estimated by the kernel, which must match the order of
derivative drv. For more information on the kernels see
section 4.3, a list of available kernels is given in A.1. The default
kernels are "MW_220" for both dimensions."model_method"
Currently available are
"iid" (independently identically distributed errors
with variance estimation from the residuals, set as default)."sarma_sep" (separable spatial ARMA (SARMA) process,
two univariate processes in both directions are estimated via
stats::arima)."sarma_HR" (fast estimation of an SARMA process by the
Hannan-Rissanen algorithm)."sarma_RSS (estimation of an separable SARMA process by
numerical minimization of the RSS)."sfarima_RSS" (estimation of a separable spatial FARIMA
(SFARIMA) model by numerical minimization of the RSS).The models and estimation methods are described in more detail in Section 4.4.
The option var_est with values
c("iid", "qarma", "sarma", "lm")has been replaced by var_model. For compatibility reasons, the option can still be used inset.optionsbut will be converted to var_model
set.options. The default values of these options typically
depend on other options and thus are put in an ellipsis. Accepted
arguments are
infl_par), the inflation exponents
(infl_exp) and trimming parameters for stabilized
estimation of the necessary derivatives (trim). Another
option to further stabilize estimation of derivatives at the boundaries
is the use of a constant estimation window at the boundaries setting the
logical flag const_window to TRUE. The default
values for the IPI-options depend partly on the regression type and
error model selected and are given below. For convenience, the contents
of IPI_options can be passed directly into the ellipsis
... of set.options().list(ar = c(1, 1), ma = c(1, 1))
(the default for SARMA and SFARIMA) specifying the model order, or any
of c("aic", "bic", "gpac") specifying an order selection
criterion. Note that gpac does not work under SFARIMA
errors.model_order. Is a list of
the form list(ar = c(1, 1), ma = c(1, 1)) (the
default).set.options() returns an object of class
"dcs_options including the following values
drv by pk = νk + 1, k = x, t.type (see default values for KR
and LP below).model_order and order_max if
available.Every argument of the
set.optionsfunction has a default value. Hence, just usingset.options()will produce a complete set of options for double conditional smoothing regression indcs(which is also implemented as default options indcs, if the argumentdcs_optionsis omitted).
Default options for kernel regression (type = "KR")
are
#> dcs_options
#> ---------------------------------------
#> options for DCS rows cols
#> ---------------------------------------
#> type: kernel regression
#> kernels used: MW_220 MW_220
#> derivative: 0 0
#> variance model:
#> ---------------------------------------
#> IPI options:
#> inflation parameters 2 1
#> inflation exponents 0.5 0.5
#> trim 0.05 0.05
#> constant window width FALSE
#> ---------------------------------------Default options for local polynomial regression
(type = "LP") are
#> dcs_options
#> ---------------------------------------
#> options for DCS rows cols
#> ---------------------------------------
#> type: local polynomial regression
#> kernel order: MW_220 MW_220
#> derivative: 0 0
#> polynomial order: 1 1
#> variance model: iid
#> ---------------------------------------
#> IPI options:
#> inflation parameters 2 1
#> inflation exponents auto
#> trim 0.05 0.05
#> constant window width FALSE
#> ---------------------------------------dcs()The dcs()-function is the main function of the package
and includes IPI-bandwidth selection and non-parametric smoothing of the
observations Y using the selected bandwidths. This function
creates an object of class dcs, which includes the results
of the DCS procedure.
Arguments of dcs() are
Y has to
have at least five rows and columns, however, for reliable results the
size should be larger."dcs_options" created by set.options(). This
argument is optional, if omitted, all options will be set to their
default values from the set.options() function."auto" if bandwidth selection should be
employed (the default).X and
T which should be ordered numerical vectors whose length
matches the number of rows of Y for X and the
number of columns of Y for T.dcs returns an object of class "dcs
including the following values
h = "auto is used in dcs, the bandwidths are
optimized via the IPI-algorithm, if h is set to fixed
values, these bandwidth are used.Y.
Either obtained by IPI bandwidth selection or given as argument in
dcs.NA, if no bandwidth selection is used.$R. Output depends on the model specification used
in dcs_options$var_model. For
var_model = "iid", it contains the estimated standard
deviation of the residuals and an indicator for stationarity, which is
true by assumption. For
dcs_options$var_model = c("sarma_sep", "sarma_RSS", "sarma_HR"),
it contains the estimated model in an object of class
"sarma" including the coefficient matrices $ar
and $ma, the standard deviation $sigma as well
as an stationarity indicator $stnry. For
dcs_options$var_model = "sfarima_RSS", the output is of
class "sfarima", with similar contents as
"sarma" and the addition of the estimated long memory
parameter vector $d."dcs_options" containing the options used in the
function.NA, if no
bandwidth selection is used.NA, if no bandwidth selection is used.surface.dcs()This function is a convenient wrapper for the
plotly::plot_ly() function of the plotly package,
for easy displaying of the considered surfaces. Direct plotting is
available for any object of class "dcs" or any numeric
matrix Y.
Arguments of surface.dcs() are
"dcs",
inheriting from a call to dcs() or a numeric matrix, which
is then directly passed to plotly::plot_ly().Y is an
object of class "dcs". Specifies the surface to be plotted,
1 for the original observations, 2 for the
smoothed surface and 3 for the residual surface. If
plot_choice is omitted and Y is an
"dcs"-object, a choice dialogue will be prompted to the
console, which asks to state one of the available options.plotly::plot_ly() function.surface.dcs() returns an object of class
"plotly".
sarma.sim()Simulation of a spatial ARMA process (SARMA). This function returns
an object of class "sarma" with attribute
"subclass" = "sim". The simulated innovations are created
from a normal distribution with specified standard deviation σ. This function uses a burn-in
period for more consistent results.
Arguments of sarma.sim() are
n_x specifies the
number of rows and n_t the number of columns. Initially, a
matrix Y' of size 2nx × 2nt
is simulated, for which simulation points with i ≤ nx
or j ≤ nt
are discarded (burn-in period).list(ar, ma, sigma). The values ar and
ma are matrices of size (px + 1) × (pt + 1)
respective (qx + 1) × (qt + 1)
and contain the coefficients in ascending lag order, so that the upper
left entry is equal to 1 (for lag 0 in both dimensions). The standard
deviation of the iid. innovations with zero mean is sigma,
which should be a single positive number. See the examples in the
application part 3.3 and the notation of SARMA models in Section
4.4.sarma.sim() returns an object of class
"sarma" with attribute "subclass" = "sim"
including the following values:
n_x, n_t).
The matrix Y is the lower left
nx × nt
submatrix of the actually simulated matrix Y′ of size 2nx × 2nt
to avoid effects from setting the initial values (burn-in period).model\$sigma). As with Y,
the original matrix has size 2nx × 2nt.sarma.est()Estimate the parameters of an SARMA of given order. It returns an
object of class "sarma" with attribute
"subclass" = "est". For estimation, three methods are
available.
Arguments of sarma.est() are
"HR"), a separable model using two univariate
estimations via stats::arima ("sep") and a
separable model which minimizes the residual sum of squares (RSS) of the
model ("RSS").list(ar = c(1, 1), ma = c(1, 1)), where all orders should
be non-negative integers. A SARMA((1, 1), (1, 1)) model is estimated by
default, if model_order is omitted.sarma.est() returns an object of class
"sarma" with attribute "subclass" = "est"
including the following values:
ar of autoregressive coefficients,
the matrix ma of moving average coefficients as well as the
standard deviation of residuals sigma.sfarima.sim()Simulation of a (separable) spatial fractional ARIMA (SFARIMA)
process. This function returns an object of class "sfarima"
with attribute "subclass" = "sim". The simulated
innovations are created from a normal distribution with specified
standard deviation σ. This
function uses a burn-in period for more consistent results.
Arguments of sfarima.sim() are
n_x specifies the
number of rows and n_t the number of columns. Initially, a
matrix Y' of size 2nx × 2nt
is simulated, for which simulation points with i ≤ nx
or j ≤ nt
are discarded (burn-in period).list(ar, ma, d, sigma).
The values ar and ma are matrices of size
(px + 1) × (pt + 1)
respective (qx + 1) × (qt + 1)
and containing the coefficients in ascending lag order, so that the
upper left entry is equal to 1 (for lag 0 in both dimensions). The
long-memory parameters dx, dt
are stored in d, a numerical vector of length 2, with 0 < dx, dt < 0.5.
The standard deviation of the iid. innovations with zero mean is
sigma, which should be a single positive number. See the
examples in the application part 3.4 and the notation of the short
memory SARMA part in Section 4.4.sfarima.sim() returns an object of class
"sfarima" with attribute "subclass" = "sim"
including the following values:
n_x, n_t).
The matrix Y is the lower left
nx × nt
submatrix of the actually simulated matrix Y′ of size 2nx × 2nt
to avoid effects from setting the initial values (burn-in period).model\$sigma). As with Y,
the original matrix has size 2nx × 2nt.sfarima.est()Estimation of an SFARIMA process. This function minimizes the
residual sum of squares (RSS) to estimate the SFARIMA-parameters of a
given order. It returns an object of class "sfarima" with
attribute "subclass" = "est".
Arguments of sfarima.est() are
list(ar = c(1, 1), ma = c(1, 1)). All orders should be
non-negative integers. A SFARIMA((1, 1), (1, 1)) model is estimated by
default, if model_order is omitted.sfarima.est() returns an object of class
"sfarima" with attribute "subclass" = "est"
including the following values:
ar of autoregressive coefficients,
the matrix ma of moving average coefficients as well as the
vector d holding the long-memory parameters and the
standard deviation of residuals sigma.kernel.assign()This function sets an external pointer to a specified boundary kernel
available in the DCSmooth package. These kernels are functions K(u, q), where
u is a vector on [q, −1] and q ∈ [0, 1]. The boundary kernels are
as proposed by Müller and Wang (1994),
Müller (1991) and constructed via the
method described in Feng and Schäfer
(2021). Available types are Müller-Wang (MW), Müller
(M) and truncated kernels (T).
Arguments of kernel.assign() are
"X_abc",
where X is specifies the type (MW,
M, T), a i the kernel order k, b is the degree of
smoothness μ and
c is the order of derivative ν. A list of currently usable kernel
identifiers can be accessed with the function
kernel.list().kernel.assign() returns an object of class
"function", which points to a precompiled kernel
function.
kernel.list()kernel.list() prints the available identifiers for use
in kernel.assign().
The argument of kernel.list() is
kernel.list() returns a list including the available
kernels as character strings, if the argument is
print = FALSE.
Available kernels are
| k | μ | ν | Truncated Kernels | Müller Kernels | Müller-Wang Kernels |
|---|---|---|---|---|---|
| 2 | 0 | 0 | M_200 |
MW_200 |
|
| 2 | 1 | 0 | M_210 |
MW_210 |
|
| 2 | 2 | 0 | T_220 |
M_220 |
MW_220 |
| 3 | 2 | 1 | T_321 |
M_321 |
MW_320 |
| 4 | 2 | 0 | T_420 |
M_420 |
MW_420 |
| 4 | 2 | 1 | M_421 |
MW_421 |
|
| 4 | 2 | 2 | T_422 |
M_422 |
MW_422 |
The DCSmooth package contains the following methods
| Function | Methods/Generics available |
|---|---|
dcs_options |
print, summary |
dcs |
plot, print,
print.summary, residuals,
summary |
sarma |
print.summary, summary |
sfarima |
print.summary, summary |
This package contains three simulated example data sets and six data sets of environmental spatial time series.
Each of the three simulated example data sets is a matrix of size 101 × 101 computed on [0, 1]2 for the following functions, where N(μ, Σ) is the bivariate normal distribution with mean vector μ and covariance matrix Σ:
The environmental application data sets feature the temperature and
wind speed surfaces of Nunn, CO (temp.nunn,
wind.nunn) and Yuma, AZ (temp.yuma,
wind.yuma). The observations are taken in 2020 in 5-minute
intervals. The temperatures are given in Celsius and wind speed
is given in m/s. All data sets consist therefore of 288 columns
(the intraday observations) and 366 rows (the days). The data is taken
from the U.S. Climate Reference Network database at www.ncdc.noaa.gov.
(see Diamond et al., 2013)
For examples of financial applications, the return and volume data of German insurance company Allianz SE is available in the package. The data is aggregated to the 5-minute level over the years 2007-2010, hence, the financial crisis 2008 is covered. The matrices consist of 1016 rows representing the days and 101 (returns) respective 102 (volumes) columns for the intraday 5-minute intervals.
The application of the package is demonstrated at the example of the
simulated function y.norm1, which represents a gaussian
peak on [0, 1]2 with nx = nt = 101
evaluation points. Different models are simulated and estimation using
dcs is demonstrated. Whenever default options are used,
they are not explicitly used as function arguments, instead only when
deviating from the defaults, the options are changed.
In order to keep the file size manageable, the
surface.dcscommands in this vignette are commented out. Run the complete code used in section 3 of this vignette bydemo("DCS_demo", package = "DCSmooth"), but note that it might take some time.
In order to set specific options use the set.options()
function to create an object of class "dcs_options".
opt1 = set.options(type = "KR", kerns = c("M_220", "M_422"), drv = c(0, 2),
var_model = "sarma_RSS",
IPI_options = list(trim = c(0.1, 0.1), infl_par = c(1, 1),
infl_exp = c(0.7, 0.7),
const_window = TRUE),
model_order = list(ar = c(1, 1), ma = c(0, 0)))
summary(opt1)
#> dcs_options
#> ---------------------------------------
#> options for DCS rows cols
#> ---------------------------------------
#> type: kernel regression
#> kernels used: M_220 M_422
#> derivative: 0 2
#> variance model:
#> ---------------------------------------
#> IPI options:
#> inflation parameters 1 1
#> inflation exponents 0.7 0.7
#> trim 0.1 0.1
#> constant window width TRUE
#> ---------------------------------------
class(opt1)
#> [1] "dcs_options"The contents of the advanced option IPI_options can be
set directly as argument in set.options(). Changing these
options might lead to non convergent bandwidths.
opt2 = set.options(type = "KR", kerns = c("M_220", "M_422"), drv = c(0, 2),
var_model = "sarma_RSS", trim = c(0.1, 0.1),
infl_par = c(1, 1), infl_exp = c(0.7, 0.7),
const_window = TRUE,
model_order = list(ar = c(1, 1), ma = c(0, 0)))
summary(opt2)
#> dcs_options
#> ---------------------------------------
#> options for DCS rows cols
#> ---------------------------------------
#> type: kernel regression
#> kernels used: M_220 M_422
#> derivative: 0 2
#> variance model:
#> ---------------------------------------
#> IPI options:
#> inflation parameters 1 1
#> inflation exponents 0.7 0.7
#> trim 0.1 0.1
#> constant window width TRUE
#> ---------------------------------------
class(opt2)
#> [1] "dcs_options"When using a model selection procedure, the additional option
order_max is available:
opt3 = set.options(var_model = "sarma_sep", model_order = "bic",
order_max = list(ar = c(0, 1), ma = c(2, 2)))
summary(opt3)
#> dcs_options
#> ---------------------------------------
#> options for DCS rows cols
#> ---------------------------------------
#> type: local polynomial regression
#> kernel order: MW_220 MW_220
#> derivative: 0 0
#> polynomial order: 1 1
#> variance model: sarma_sep
#> ---------------------------------------
#> IPI options:
#> inflation parameters 2 1
#> inflation exponents auto
#> trim 0.05 0.05
#> constant window width FALSE
#> ---------------------------------------The example data set is simulated by using iid. errors:
# surface.dcs(y.norm1)
y_iid = y.norm1 + matrix(rnorm(101^2), nrow = 101, ncol = 101)
# surface.dcs(y_iid)While local linear regression has some clear advantages over kernel regression, kernel regression is the faster method.
opt_iid_KR = set.options(type = "KR")
dcs_iid_KR = dcs(y_iid, opt_iid_KR)
# print results
dcs_iid_KR
#> dcs
#> ------------------------------------------
#> DCS with automatic bandwidth selection:
#> ------------------------------------------
#> Selected Bandwidths:
#> h_x: 0.18875
#> h_t: 0.19282
#> Variance Factor:
#> c_f: 1
#> ------------------------------------------
# print options used for DCS procedure
dcs_iid_KR$dcs_options
#> dcs_options
#> ---------------------------------------
#> options for DCS rows cols
#> ---------------------------------------
#> type: kernel regression
#> kernels used: MW_220 MW_220
#> derivative: 0 0
#> variance model: iid
#> ---------------------------------------
# plot regression surface
# surface.dcs(dcs_iid_KR, plot_choice = 2)The summary of the "dcs"-object provides some more
detailed information:
summary(dcs_iid_KR)
#> summary_dcs
#> ------------------------------------------
#> DCS with automatic bandwidth selection:
#> ------------------------------------------
#> Results of kernel regression:
#> Estimated Bandwidths: h_x: 0.1888
#> h_t: 0.1928
#> Variance Factor: c_f: 1
#> Iterations: 4
#> Time used (seconds): 0.009931
#> ------------------------------------------
#> Variance Model: iid
#> ------------------------------------------
#> sigma: 0.9968943
#> stationary: TRUE
#> ------------------------------------------
#> See used parameter with "$dcs_options".This is the default method, specification of options is not necessary. Note that local polynomial regression requires the bandwidth to cover at least the number of observations of the polynomial order plus one. For small bandwidths or too few observation points in one dimension, local polynomial regression might fail (“Bandwidth h must be larger for local polynomial regression.”). It is suggested to use kernel regression in this case.
dcs_LP_iid = dcs(y_iid)
dcs_LP_iid
#> dcs
#> ------------------------------------------
#> DCS with automatic bandwidth selection:
#> ------------------------------------------
#> Selected Bandwidths:
#> h_x: 0.1941
#> h_t: 0.20778
#> Variance Factor:
#> c_f: 1.0003
#> ------------------------------------------
summary(dcs_LP_iid)
#> summary_dcs
#> ------------------------------------------
#> DCS with automatic bandwidth selection:
#> ------------------------------------------
#> Results of local polynomial regression:
#> Estimated Bandwidths: h_x: 0.1941
#> h_t: 0.2078
#> Variance Factor: c_f: 1
#> Iterations: 4
#> Time used (seconds): 0.02297
#> ------------------------------------------
#> Variance Model: iid
#> ------------------------------------------
#> sigma: 0.9970611
#> stationary: TRUE
#> ------------------------------------------
#> See used parameter with "$dcs_options".
# plot regression surface
# surface.dcs(dcs_LP_iid, plot_choice = 2)A matrix containing innovations following a SARMA((px, pt), (qx, qt))
process can be obtained by using the sarma.sim() function.
We use the following SARMA((1, 1), (1, 1))-process as example: $$
AR = \begin{pmatrix} 1 & 0.4 \\ -0.3 & -0.12 \end{pmatrix},
\quad
MA = \begin{pmatrix} 1 & -0.5 \\ -0.2 & 0.1 \end{pmatrix}
\quad \text{ and } \quad
\sigma^2 = 0.25
$$
ar_mat = matrix(c(1, -0.3, 0.4, 0.12), nrow = 2, ncol = 2)
ma_mat = matrix(c(1, -0.2, -0.5, 0.1), nrow = 2, ncol = 2)
sigma = sqrt(0.25)
model_list = list(ar = ar_mat, ma = ma_mat, sigma = sigma)
sim_sarma = sarma.sim(n_x = 101, n_t = 101, model = model_list)
# SARMA observations
y_sarma = y.norm1 + sim_sarma$Y
# surface.dcs(y_sarma)Estimation of an SARMA process for a given order is implemented via
the sarma.est() function (note that the simulated matrix
can be accessed via $Y):
est_sarma = sarma.est(sim_sarma$Y, method = "HR",
model_order = list(ar = c(1, 1), ma = c(1, 1)))
summary(est_sarma)
#> ------------------------------------------
#> Estimation of SARMA((1,1),(1,1))
#> ------------------------------------------
#> sigma: 0.5008
#> stationary: TRUE
#> ar:
#> lag 0 lag 1
#> lag 0 1.0000 0.4042
#> lag 1 -0.3051 0.1291
#>
#> ma:
#> lag 0 lag 1
#> lag 0 1.0000 -0.49690
#> lag 1 -0.1804 0.07785The dcs()-command is used with the default SARMA((1, 1), (1, 1)) model (correctly specified)
and with an SARMA((1, 1), (0, 0))
(i.e. SAR(1, 1)) model. The chosen
estimation procedure is "sep":
# SARMA((1, 1), (1, 1))
opt_sarma_1 = set.options(var_model = "sarma_sep")
dcs_sarma_1 = dcs(y_sarma, opt_sarma_1)
summary(dcs_sarma_1$var_est)
#> ------------------------------------------
#> Estimation of SARMA((1,1),(1,1))
#> ------------------------------------------
#> sigma: 0.5441
#> stationary: TRUE
#> ar:
#> lag 0 lag 1
#> lag 0 1.0000 0.5003
#> lag 1 -0.5246 -0.2624
#>
#> ma:
#> lag 0 lag 1
#> lag 0 1.0000 -0.47000
#> lag 1 -0.1524 0.07165
# SARMA((1, 1), (0, 0))
# Use "sarma_HR" as it includes fast Yule-Walker estimation
opt_sarma_2 = set.options(var_model = "sarma_HR",
model_order = list(ar = c(1, 1), ma = c(0, 0)))
dcs_sarma_2 = dcs(y_sarma, opt_sarma_2)
summary(dcs_sarma_2$var_est)
#> ------------------------------------------
#> Estimation of SARMA((1,1),(0,0))
#> ------------------------------------------
#> sigma: 0.5331
#> stationary: TRUE
#> ar:
#> lag 0 lag 1
#> lag 0 1.000 0.63580
#> lag 1 -0.176 0.09775
#>
#> ma:
#> lag 0
#> lag 0 1Automated order selection is used with
model_order = c("aic", "bic", "gpac") in
set.options(). The first one minimizes the AIC, the second
one the BIC and the third uses a generalized partial autocorrelation
function for order selection (not available for SFARIMA estimation).
Order selection for large data sets is slowly in general, however, the
"gpac" might be slightly faster than the other two. If
automatic order selection is used, the argument order_max
sets the maximum orders to be tested in the same way as
model_order is usually defined as a list.
# BIC
opt_sarma_3 = set.options(var_model = "sarma_HR", model_order = "bic",
order_max = list(ar = c(2, 2), ma = c(2, 2)))
dcs_sarma_3 = dcs(y_sarma, opt_sarma_3)
summary(dcs_sarma_3$var_est)
#> ------------------------------------------
#> Estimation of SARMA((1,1),(1,1))
#> ------------------------------------------
#> sigma: 0.4997
#> stationary: TRUE
#> ar:
#> lag 0 lag 1
#> lag 0 1.0000 0.4031
#> lag 1 -0.3083 0.1271
#>
#> ma:
#> lag 0 lag 1
#> lag 0 1.0000 -0.50190
#> lag 1 -0.1879 0.07439
# gpac
opt_sarma_4 = set.options(var_model = "sarma_HR", model_order = "gpac",
order_max = list(ar = c(2, 2), ma = c(2, 2)))
dcs_sarma_4 = dcs(y_sarma, opt_sarma_4)
summary(dcs_sarma_4$var_est)
#> ------------------------------------------
#> Estimation of SARMA((2,2),(2,0))
#> ------------------------------------------
#> sigma: 0.5061
#> stationary: TRUE
#> ar:
#> lag 0 lag 1 lag 2
#> lag 0 1.00000 0.85200 0.291300
#> lag 1 -0.46200 -0.14640 -0.051140
#> lag 2 0.01744 -0.02795 -0.009702
#>
#> ma:
#> lag 0
#> lag 0 1.00000
#> lag 1 -0.33880
#> lag 2 -0.01287This package includes a bandwidth selection algorithm when the errors ε(x, t) follow an SFARIMA((px, pt), (qx, qt)) process with long memory. Order selection for SFARIMA models works exactly as in the SARMA case. We use the same SARMA model as in 3.3 with long-memory parameters d = (0.3, 0.1):
ar_mat = matrix(c(1, -0.3, 0.4, 0.12), nrow = 2, ncol = 2)
ma_mat = matrix(c(1, -0.2, -0.5, 0.1), nrow = 2, ncol = 2)
d = c(0.3, 0.1)
sigma = sqrt(0.25)
model_list = list(ar = ar_mat, ma = ma_mat, d = d, sigma = sigma)
sim_sfarima = sfarima.sim(n_x = 101, n_t = 101, model = model_list)
# SFARIMA surface observations
y_sfarima = y.norm1 + sim_sfarima$Y
# surface.dcs(y_sfarima)
opt_sfarima = set.options(var_model = "sfarima_RSS")
dcs_sfarima = dcs(y_sfarima, opt_sfarima)
summary(dcs_sfarima$var_est)
#> ------------------------------------------
#> Estimation of SFARIMA((1,1),(1,1))
#> ------------------------------------------
#> d: 0.3146 0.09794
#> SD (sigma): 0.4978
#> stationary: TRUE
#> ar:
#> lag 0 lag 1
#> lag 0 1.0000 0.38940
#> lag 1 -0.1575 -0.06132
#>
#> ma:
#> lag 0 lag 1
#> lag 0 1.00000 -0.50950
#> lag 1 -0.07279 0.03708Local polynomial estimation is suitable for estimation of derivatives
of a function or a surface. While estimation of derivatives works as
well under dependent errors, the example uses the iid. model from 3.2.
Derivatives can be computed for any derivative vector drv,
if the values are non-negative and an appropriate kernel function is
chosen (such that the derivative order of the kernel matches the
derivative order in drv). Note that the order of the
polynomials for the νth
derivative is chosen to be pi = νi + 1, i = x, t.
As bandwidths increase with the order of the derivatives, the bandwidth
might be large for higher derivative orders.
The estimator for m(1, 0)(x, t) is
opt_drv_1 = set.options(drv = c(1, 0), kerns = c("MW_321", "MW_220"))
opt_drv_1$IPI_options$trim = c(0.1, 0.1)
dcs_drv_1 = dcs(y_iid, opt_drv_1)
dcs_drv_1
#> dcs
#> ------------------------------------------
#> DCS with automatic bandwidth selection:
#> ------------------------------------------
#> Selected Bandwidths:
#> h_x: 0.24868
#> h_t: 0.41024
#> Variance Factor:
#> c_f: 1.002
#> ------------------------------------------
# surface.dcs(dcs_drv_1, trim = c(0.1, 0.1), plot_choice = 2)The estimator for m(0, 2)(x, t) is given by
opt_drv_2 = set.options(drv = c(0, 2), kerns = c("MW_220", "MW_422"))
opt_drv_2$IPI_options$trim = c(0.1, 0.1)
dcs_drv_2 = dcs(y_iid, opt_drv_2)
dcs_drv_2
#> dcs
#> ------------------------------------------
#> DCS with automatic bandwidth selection:
#> ------------------------------------------
#> Selected Bandwidths:
#> h_x: 0.28678
#> h_t: 0.37157
#> Variance Factor:
#> c_f: 1.002
#> ------------------------------------------
# surface.dcs(dcs_drv_2, trim = c(0.1, 0.1), plot_choice = 2)The double conditional smoothing (DCS, see Feng (2013)) is a spatial smoothing technique which effectively reduces the two-dimensional estimation to two one-dimensional estimation procedures. The DCS is defined for kernel regression as well as for local polynomial regression.
Classical bivariate (and multivariate) regression has been considered e.g. by Herrmann et al. (1995) (kernel regression) and Ruppert and Wand (1994) (local polynomial regression). The DCS provides now a faster and, especially for equidistant data, more efficient smoothing scheme, which leads to reduced computation time. For the DCS procedure implemented in this package, consider a (nx × nt)-matrix Y of non-empty observations ui, j and equidistant covariates X, T on [0, 1], where X has length nx and T has length nt. The model is then yi, j = m(xi, tj) + εi, j where m(x, t) is the mean or trend function, xi ∈ X, tj ∈ T and ε is a random error function with zero mean. The model in matrix form is Y = M1 + E at the observation points.
The main assumption of the DCS is that of product kernels, i.e. the weights in the respective methods are constructed by K(u, v) = K1(u)K2(v). Now, a two stage smoother can be constructed by either the kernel weights directly (kernel regression) or by using locally weighted regression with kernels K1, K2, in any case, the weights are called Wx and Wt. The DCS procedure implemented in DCSmooth smoothes over rows (conditioning on X) first and then over columns (conditioning on T), although switching the smoothing order is exactly equivalent. Hence, the DCS is given by the following equations: $$ \begin{align*} \mathbf{\widehat{M}}_0[ \ , j] &= \mathbf{Y} \cdot \mathbf{W}_t[ \ , j] \\ \mathbf{\widehat{M}}_1[i, \ ] &= \mathbf{W}_x [i, \ ] \cdot \mathbf{\widehat{M}}_0 \end{align*} $$
The bandwidth vector h = (hx, ht) is selected via an iterative plug-in (IPI) algorithm (Gasser et al., 1991). The IPI selects the optimal bandwidths by minimizing the mean integrated squared error (MISE) of the estimator. As the MISE includes derivatives of the regression surface m(x, t), auxiliary bandwidths for estimation of these derivatives are calculated via an inflation method. These inflation method connects the bandwidths of m(x, t) with that of a derivative m(νx, νt)(x, t) by h̃k = ck ⋅ hkα, k = x, t and is called exponential inflation method (EIM). The values of ck are chosen on simulations, that of alpha are subject to the derivative of interest. The IPI now starts with an initial bandwidth h0 (chosen to be h0 = (0.1, 0.1)) and calculates in each step s the auxiliary bandwidths h̃k, s from hs − 1 and hs from the smoothed derivative surfaces using h̃k, s. The iteration process finishes until a certain threshold is reached.
In kernel regression, the boundary problem exists, which leads to biased estimated at the boundaries of the regression surface. This problem can (partially) be solved by means of suitable boundary kernels as introduced by Müller (1991) and Müller and Wang (1994). These boundary kernels differ in their degrees of smoothness and hence lead to different estimation results at the boundaries. However, all kernels are similar to the classical kernels in the interior region of the regression.
Following Feng and Schäfer (2021), a
boundary modification is also defined for local polynomial regression.
In the DCSmooth package, the local polynomial regression is
always with boundary modification weights. Kernel types available
(either for kernel regression or local polynomial regression) are
Müller-type, Müller-Wang-type and truncated kernels, denoted by
M, MW and T. In most
applications, the Müller-Wang type are the preferred weighting
functions.
For observations X = xi,
xi, i = 1, …, nx
and a given bandwidth h,
define $u_r = \frac{x_r - X}{h} \in [-1,
q]$. A (left) boundary kernel function Kql(u)
of order (k, μ, ν) is
defined on [−1, q], for q ∈ [0, 1] and has the following
properties $$
\int_{-1}^q u^j K^l_q(u) \text{d} u = \begin{cases}
0 & \text{ for } 0 \leq j < k, j \neq \nu\\
(-1)^\nu \nu! & \text{ for } j = \nu\\
\beta \neq 0 & \text{ for } j = k
\end{cases}
$$ The corresponding right boundary kernels can be calculated by
Kqr(u) = (−1)νKq(−u).
The boundary kernels assigned by kernel.assign() are left
boundary kernels.
The SARMA process εi, j
is given by the following equations: $$
\phi(B_1, B_2)\varepsilon_{i,j} = \psi(B_1, B_2)\eta_{i,j},\\
$$ where the lag operators are B1εi, j = εi − 1, j
and B2εi, j = εi, j − 1,
$\xi \underset{iid.}{\sim}
\mathcal{N}(0,\sigma^2)$ and $$
\phi(z_1, z_2) = \sum_{m = 0}^{p_1} \sum_{n = 0}^{p_2} \phi_{m, n} z_1^m
z_2^n, \ \psi(z_1, z_2) = \sum_{m = 0}^{q_1} \sum_{n = 0}^{q_2} \psi_{m,
n} z_1^m z_2^n.
$$ The coefficients ψm, n
and ϕm, n
are written in matrix form $$
\begin{align*}
\boldsymbol{\phi} = \begin{pmatrix}
\phi_{0,0} & \dots & \phi_{0, p_2} \\
\vdots & \ddots & \\
\phi_{p_1, 0} & & \phi_{p_1, p_2}
\end{pmatrix} \ \text{ and } \ \boldsymbol{\psi} = \begin{pmatrix}
\psi_{0, 0} & \dots & \psi_{0, q_2} \\
\vdots & \ddots & \\
\psi_{q_1, 0} & & \psi_{q_1, q_2}
\end{pmatrix},
\end{align*}
$$ where Φ is the
AR-part ($var_model$ar) and Ψ is the MA-part
($var_model$ma). The example from 3.3, $$
\begin{align}
\boldsymbol{\phi} = \begin{pmatrix} 1 & 0.4 \\ -0.3 & -0.12
\end{pmatrix}
\ \text{ and } \
\boldsymbol{\psi} = \begin{pmatrix} 1 & -0.2 \\ -0.5 & 0.1
\end{pmatrix}
\end{align},
$$ would then reduce to the process εi, j = 0.4εi, j − 1 − 0.3εi−, j + 0.2εi − 1, j − 1 + 0.2ξi, j − 1 + 0.2ξi − 1, j − 0.5ξi − 1, j − 1 + ξi, j.
Note that this process can be written as product of two univariate
processes in the sense that $$
\phi_1(B_1)\phi_2(B_2)^T \varepsilon_{i,j} = \psi_1(B_1) \psi_2(B_2)^T
\eta_{i,j},\\
$$ with $$
\phi_1 = \begin{pmatrix} 1 \\ -0.3 \end{pmatrix}, \quad \phi_2 =
\begin{pmatrix} 1 \\ 0.4 \end{pmatrix}, \quad \psi_1 = \begin{pmatrix} 1
\\ -0.5 \end{pmatrix}, \quad \psi_2 = \begin{pmatrix} 1 \\ -0.2
\end{pmatrix}.
$$ Hence, these process forms a separable SARMA. Estimation of
separable SARMA models can be reduced to estimation of univariate ARMA
models.
For estimation of SARMA models, three methods are implemented in DCSmooth:
This method is only available under the assumption of a separable model. Define two univariate time series ε1, r = {ε1}r = {εi, j}i + nt(j − 1), ε2, s = {ε2}s = {εi, j}j + nx * (i − 1) for r, s = 1, …, nx ⋅ nt, i = 1, …, nx, j = 1, …, nt. The parameters ϕ1, ψ1 of ε1, r and ϕ2, ψ2 of ε2, s can then be estimated by well-known maximum likelihood estimators.
The SARMA model can be rewritten as $$ \eta_{i,j} = \psi(B_1, B_2)^{-1}\phi(B_1, B_2)\varepsilon_{i,j},\\ $$ which allows for an AR(∞)-representation of the SARMA model $$ \eta_{i,j} = \sum_{r,s = 0}^{\infty} \theta^{AR}_{r,s} \varepsilon_{i-r, s-j} $$
From this, we can define the residual sum of squares (RSS) and get an estimate for the vector of coefficients θ = c(ϕ1, 0, ϕ0, 1, …, ψ1, 0, ψ0, 1, … by $$ \widehat{\theta} = \arg \min RSS \approx \arg \min \sum_{i,j = 0}^{\infty} \eta_{i,j}^2 $$
Calculation of the $AR() representation of an SARMA model is difficult for a general SARMA but for a separable SARMA, the known univariate formulas hold.
These procedure can be directly used for SFARIMA models, if the long memory parameter d is included in θ (see Beran et al. (2009)).
The previously defined estimation methods require numeric optimization of some quantities and hence take more time for calculation on a computer. The Hannan-Rissanen algorithm (Hannan and Rissanen (1982)) provides a much faster estimation procedure. An extension to SARMA models is straightforward:
The main idea Hannan-Rissanen algorithm is to use a high-order SAR auxiliary model for initial estimation of the non observable innovations sequence. Then, a linear regression model is applied, which yields the SARMA-parameters from the observations and estimated innovations.
Let {ε}i, j be the the ordered observations of the SARMA((px, pt), (qx, qt)) process ε(x, t) with i = 1, …, nx, j = 1, …, nt. The SARMA parameters $\mathbf{\phi, \mathbf{\psi}$ are then estimated by the modified Hannan-Rissanen algorithm: 1. Obtain the auxiliary residuals η̃i, j by fitting a high-order autoregressive model with (p̃x, p̃t) ≥ (px, pt) to the observations: $$ \begin{align} \widetilde{\eta}_{i,j} = \begin{cases} \sum_{m = 0}^{\widetilde{p}_1} \sum_{n = 0}^{\widetilde{p}_2} \widetilde{\phi}_{m,n} \varepsilon_{i - m, j - n}, & \widetilde{p}_1 < i \leq n_1, \widetilde{p}_2 < j \leq n_2\\ 0, &1 \leq i \leq \widetilde{p}_1, 1 \leq j \leq \widetilde{p}_2 \end{cases} \end{align} $$ where ϕ̃m, n is estimated by the Yule-Walker equations and ϕ0, 0 = 1. 2. Obtain ϕ̂m, n and ψ̂m, n and the estimated innovations η̂i, j by linear regression from $$ \begin{align} \varepsilon_{i,j} = -\sum_{\substack{m=0\\m \neq n = 0}}^{p_1} \sum_{\substack{n=0\\n \neq m = 0}}^{p_2} \widehat{\phi}_{m,n} \varepsilon_{i-m,j-n} + \sum_{\substack{m=0\\m \neq n = 0}}^{q_1} \sum_{\substack{n=0\\n \neq m = 0}}^{q_2} \widehat{\psi}_{m,n} \widetilde{\xi}_{i-m,j-n} + \widehat{\xi}_{i,j} \end{align} $$ The resulting coefficients $\widehat{\mathbf{\phi}}$, $\widehat{\mathbf{\psi}}$ are then the estimates for the parameters.
The autocovariance function of the SARMA-process is γ(s, t) = 𝔼(εi, jεi + s, j + t). For p̃x, p̃t, the spatial Yule-Walker equation for the SAR(p̃x, p̃t) approximation of the SARMA is then $$ \begin{align} \mathbf{\Gamma} \ \left( \text{vec} (\widetilde{\boldsymbol{\phi}}) \setminus \{ \widetilde{\phi}_{0,0}\} \right) = \gamma \end{align} $$ where Γ is the autocovariance matrix and γ the vector of autocovariances when using the following notation: $$ \begin{align} \mathbf{\Gamma} &= \begin{pmatrix} \gamma(0,0) & \gamma(1,0) & \dots & \gamma(\widetilde{p}_1,0) & \gamma(2,0) & \dots & \gamma(\widetilde{p}_1 - 1, \widetilde{p}_2)\\ \gamma(1,0) & \gamma(0,0) & & & & & \gamma(\widetilde{p}_1 - 2, \widetilde{p}_2)\\ \vdots & & \ddots & & & & \vdots \\ \gamma(\widetilde{p}_1 - 1, \widetilde{p}_2) & & \dots & & & & \gamma(0,0) \end{pmatrix} \\ \gamma &= \left(\gamma(0,1), \gamma(0,2), \dots, \gamma(\widetilde{p}_1, 0), \dots, \gamma(\widetilde{p}_1, \widetilde{p}_2) \right)^T. \end{align} $$