.mom
file as a gnssts
objectLet us first load the gmwmx
package.
library(gmwmx)
dobs = gmwmx::PBO_get_station("DOBS", column = "dN") write.gnssts(dobs, filename = "data_dobs.mom")
Consider that you want to estimate a model on data saved in a .mom
file located at a specific file_path
on your computer, where file_path
is the path where is located the .mom
file (for example file_path = "/home/name_of_the_user/Documents/data.mom"
)
file_path = system.file("extdata", "data_dobs.mom", package = "gmwmx", mustWork = T) data_dobs = read.gnssts(filename = file_path)
For example, the corresponding .mom
file could have a similar looking:
# sampling period 1.000000 # offset 55285.000000 # offset 58287.770833 52759.5 -0.01165 52760.5 -0.01102 52761.5 -0.01147 ...
You can import the .mom
file as a with the function read.gnssts()
as such:
data_dobs = read.gnssts(filename = file_path)
Objects created or imported with create.gnss()
or read.gnssts()
are of class gnssts
.
class(data_dobs)
By inspecting the structure of a gnssts
object, we observe that gnssts
objects specify the time vector, the observation vector, the sampling period and the times at which there are location shifts (jumps).
str(data_dobs)
We can represent the signal as such:
plot(data_dobs$t, data_dobs$y, type="l")
The gmwmx
package allows to estimate linear model with correlated residuals that are described by a functional model and a stochastic noise model.
More precisely, for the functional model, we consider a linear model which can be expressed as:
\begin{equation} \mathbf{Y} = \mathbf{A} {{\bf x}}_0 + \boldsymbol{\varepsilon},
\end{equation}
where $\mathbf{Y} \in {\rm I!R}^n$ denotes the response variable of interest (i.e., vector of GNSS observations), $\mathbf{A} \in {\rm I!R}^{n \times p}$ a fixed design matrix, ${{\bf x}}_0 \in \mathcal{X} \subset {\rm I!R}^p$ a vector of unknown constants and $\boldsymbol{\varepsilon} \in {\rm I!R}^n$ a vector of (zero mean) residuals.
The gmwmx
package allows to estimate functional models for which the $i$-th component of the vector $\mathbf{A} {{\bf x}}_0$ can be described as follows:
\begin{equation} \mathbb{E}[\mathbf{Y}i] = \mathbf{A}_i^T {{\bf x}}_0 = a+b\left(t{i}-t_{0}\right)+\sum_{h=1}^{2}\left[c_{h} \sin \left(2 \pi f_{h} t_{i}\right)+d_{h} \cos \left(2 \pi f_{h} t_{i}\right)\right] + \sum_{k=1}^{n_{g}} g_{k} H\left(t_{i}-t_{k}\right), \end{equation}
where $a$ is the initial position at the reference epoch $t_0$, $b$ is the velocity parameter, $c_k$ and $d_k$ are the periodic motion parameters ($h = 1$ and $h = 2$ represent the annual and semi-annual seasonal terms, respectively). The offset terms models earthquakes, equipment changes or human intervention in which $g_k$ is the magnitude of the change at epochs $t_k$, $n_g$ is the total number of offsets, and $H$ is the Heaviside step function. Note that the estimates of the parameters of the functional model are provided in unit/day.
Regarding the stochastic model, we assume that $\boldsymbol{\varepsilon}_i=\mathbf{Y}_i-\mathbb{E}[\mathbf{Y}_i]$ is a strictly (intrinsically) stationary process and that
\begin{equation} \boldsymbol{\varepsilon} \sim \mathcal{F} \left{\mathbf{0}, \boldsymbol{\Sigma}(\boldsymbol{\gamma}_0)\right} , \label{eq:model:noise} \end{equation}
where $\mathcal{F}$ denotes some probability distribution in ${\rm I!R}^n$ with mean ${\bf 0}$ and covariance $\boldsymbol{\Sigma}(\boldsymbol{\gamma}_0)$.
We assume that $\boldsymbol{\Sigma}(\boldsymbol{\gamma}_0) > 0$ and that it depends on the unknown parameter vector $\boldsymbol{\gamma}_0 \in \boldsymbol{\Gamma} \subset {\rm I!R}^q$. This parameter vector specifies the covariance of the observations and is often referred to as the stochastic parameters.
Hence, we let $\boldsymbol{\theta}_0 = \left[\boldsymbol{{\bf x}}_0^{\rm T} \;\; \boldsymbol{\gamma}_0^{\rm T}\right]^{\rm T} \in \boldsymbol{\Theta} = \mathcal{X} \times \boldsymbol{\Gamma} \subset {\rm I!R}^{p + k}$ denote the unknown parameter vector of the model described above.
The gmwmx
allows to estimate parameters of a specified functional model as well as parameters of a stochastic model (i.e. $\hat{\boldsymbol{\theta}} = \left[\boldsymbol{\hat{\boldsymbol{x}}}^{T} \;\; \hat{\boldsymbol{\gamma}}^{T}\right]$) defined by a combinations of
(1) White noise (wn
)
(2) Matérn process (matern
)
(3) Fractional Gaussian noise (fgn
) and
(4) Power Law process (powerlaw
).
Note that only the gmwmx
current version accepts only one process of each kind.
You can estimate a model using the GMWMX estimator with the function estimate_gmwmx()
.
The stochastic model considered is specified by a string provided to the argument model_string
which is a combination of the strings wn
, powerlaw
, matern
and fgn
separated by the character +
.
You specify the initialization values for solving the optimization problem at the GMWM estimation step that estimate the stochastic model by providing a numeric vector of the correct length (the total number of parameters of the stochastic model specified in model_string
) to the argument theta_0
.
You can compute confidence intervals for estimated functional parameters of an estimated model by setting the argument ci
to TRUE
.
Let us consider a single sinusoidal signal with the jumps specified in the gnssts
object and a combination of a White noise and a Power Law process for the stochastic model.
fit_dobs_wn_plp_gmwmx = estimate_gmwmx(x = data_dobs, theta_0 = c(0.1, 0.1, 0.1), model_string = "wn+powerlaw", n_seasonal = 1, ci = T)
file_path = system.file("extdata", "fit_dobs_wn_plp_gmwmx.rda", package = "gmwmx", mustWork = T) load(file_path)
Estimated models are of class gnsstsmodel
class(fit_dobs_wn_plp_gmwmx)
We can print the estimated model or extract estimated parameters (functional and stochastic) as such:
print(fit_dobs_wn_plp_gmwmx) fit_dobs_wn_plp_gmwmx$beta_hat fit_dobs_wn_plp_gmwmx$theta_hat
We can also plot graphically the estimated functional model on the time series and the Wavelet variance of residuals by calling the plot.gnsstsmodel
method on a gnsstsmodel
object.
plot(fit_dobs_wn_plp_gmwmx)
We can specify the number of iterations of the GMWMX to compute respectively the GMWMX-1 and GMWMX-2 or other iteration of the GMWMX with the argument k_iter
. For example we can compute the GMWMX-2 as such:
fit_dobs_wn_plp_gmwmx_2 = estimate_gmwmx(x = data_dobs, theta_0 = c(0.1, 0.1, 0.1), model_string = "wn+powerlaw", n_seasonal = 1, k_iter = 2)
Assuming that you have Hector available on the PATH
, an estimation of the model can the be performed using the Maximum Likelihood Estimation (MLE) method implemented in Hector as such:
fit_dobs_wn_plp_mle = estimate_hector(x = data_dobs, model_string = "wn+powerlaw", n_seasonal = 1)
file_path_mle = system.file("extdata", "fit_dobs_wn_plp_mle.rda", package = "gmwmx", mustWork = T) load(file_path_mle)
Similarly we can plot and extract the model parameters of the estimated model:
plot(fit_dobs_wn_plp_mle) fit_dobs_wn_plp_mle$beta_hat fit_dobs_wn_plp_mle$theta_hat
We can load time series data from the Plate Boundary Observatory (PBO) as gnssts
object with PBO_get_station()
:
cola = PBO_get_station("COLA", column = "dE")
save(cola, file="cola.rda")
cola_path = system.file("extdata", "cola.rda", package = "gmwmx", mustWork = T) load(cola_path)
Let us consider three potential models for the stochastic model of this signal. More precisely let us consider:
fit_cola_wn_plp = estimate_gmwmx(cola, model_string = "wn+powerlaw", theta_0 = c(0.1,0.1,0.1), n_seasonal = 1, ci = T)
file_path = system.file("extdata", "fit_cola_wn_plp.rda", package = "gmwmx", mustWork = T) load(file_path)
plot(fit_cola_wn_plp)
fit_cola_wn_fgn = estimate_gmwmx(cola, model_string = "wn+fgn", theta_0 = c(0.1,0.1,0.2), n_seasonal = 1, ci = T)
file_path = system.file("extdata", "fit_cola_wn_fgn.rda", package = "gmwmx", mustWork = T) load(file_path)
plot(fit_cola_wn_fgn)
fit_cola_wn_matern = estimate_gmwmx(cola, model_string = "wn+matern", theta_0 = c(0.1,0.1,0.1,0.1), n_seasonal = 1, ci = T)
file_path = system.file("extdata", "fit_cola_wn_matern.rda", package = "gmwmx", mustWork = T) load(file_path)
plot(fit_cola_wn_matern)
You can compare estimated models with the function compare_fits()
compare_fits(fit_cola_wn_plp, fit_cola_wn_matern)
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.