Package 'gmwmx'

Title:	Estimate Functional and Stochastic Parameters of Linear Models with Correlated Residuals
Description:	Implements the Generalized Method of Wavelet Moments with Exogenous Inputs estimator (GMWMX) presented in Cucci, D. A., Voirol, L., Kermarrec, G., Montillet, J. P., and Guerrier, S. (2023) <doi:10.1007/s00190-023-01702-8>. The GMWMX estimator allows to estimate functional and stochastic parameters of linear models with correlated residuals. The 'gmwmx' package provides functions to estimate, compare and analyze models, utilities to load and work with Global Navigation Satellite System (GNSS) data as well as methods to compare results with the Maximum Likelihood Estimator (MLE) implemented in Hector.
Authors:	Davide Antonio Cucci [aut], Lionel Voirol [aut, cre], Stéphane Guerrier [aut], Jean-Philippe Montillet [ctb], Gaël Kermarrec [ctb]
Maintainer:	Lionel Voirol <[email protected]>
License:	AGPL-3
Version:	1.0.3
Built:	2025-02-05 04:41:20 UTC
Source:	https://github.com/cran/gmwmx

Help Index

GNSS time series for PBO station COLA
Compare graphically two gnsstsmodel objects.
Define matrix A of the functional model
Create a gnssts object
Estimate a stochastic model in a two-steps procedure using the GMWMX estimator.
Estimate a stochastic model based on the MLE and the Hector implementation.
Extract offsets for a PBO station
Load station data from PBO
Plotting method for a gnsstsmodel object.
Print method for a gnsstsmodel object.
Read a gnssts object
Remove outliers from a gnssts object using Hector
Write a gnssts object

GNSS time series for PBO station COLA

Description

Data from station COLA of the Plate Boundary Observatory

Usage

cola
cola

Format

A gnssts object of the East position (dE) of the COLA station

Source

https://data.unavco.org/archive/gnss/products/position/COLA/COLA.pbo.igs14.pos

Compare graphically two `gnsstsmodel` objects.

Description

Compare graphically two gnsstsmodel objects.

Usage

compare_fits(fit_1, fit_2, main = NULL, y_unit = "mm", x_unit = "days")
compare_fits(fit_1, fit_2, main = NULL, y_unit = "mm", x_unit = "days")

Arguments

`fit_1`	A `gnsstsmodel` object.
`fit_2`	A `gnsstsmodel` object.
`main`	A `string` specifying the plot title.
`y_unit`	A `string` specifying the y axis label.
`x_unit`	A `string` specifying the x axis label.

Value

No return value. Produce a plot comparing two estimated models.

Examples

## Not run: 
data(cola)
fit_gmwmx_1 = estimate_gmwmx(x = cola,
                             theta_0 = c(0.1,0.1,0.1,0.1), 
                             n_seasonal = 1, 
                             model_string = "wn+matern")
fit_gmwmx_2 = estimate_gmwmx(x = cola,
                             theta_0 = c(0.1,0.1,0.1), 
                             n_seasonal = 1, 
                             model_string = "wn+powerlaw")
compare_fits(fit_gmwmx_1, fit_gmwmx_2)

## End(Not run)
## Not run: 
data(cola)
fit_gmwmx_1 = estimate_gmwmx(x = cola,
                             theta_0 = c(0.1,0.1,0.1,0.1), 
                             n_seasonal = 1, 
                             model_string = "wn+matern")
fit_gmwmx_2 = estimate_gmwmx(x = cola,
                             theta_0 = c(0.1,0.1,0.1), 
                             n_seasonal = 1, 
                             model_string = "wn+powerlaw")
compare_fits(fit_gmwmx_1, fit_gmwmx_2)

## End(Not run)

Define matrix A of the functional model

Description

Define matrix A of the functional model

Usage

create_A_matrix(t_nogap, jumps, n_seasonal)
create_A_matrix(t_nogap, jumps, n_seasonal)

Arguments

`t_nogap`	A `vector` specifying the index of the time series.
`jumps`	A `vector` specifying the time at which there is a mean shift of the time series. Should be specified to `NULL` if there is not presence of offsets in the signal.
`n_seasonal`	An `integer` specifying the number of sinusoidal signals in the time series.

Value

Matrix A in order to compute the functional component of the model in a linear fashion

Examples

n= 10*365
jump_vec <- c(200, 300, 500)
nbr_sin = 2
A <- create_A_matrix(1:n, jump_vec, n_seasonal = nbr_sin)
head(A)
A <- create_A_matrix(1:n, jumps = NULL, n_seasonal = nbr_sin)
head(A)
n= 10*365
jump_vec <- c(200, 300, 500)
nbr_sin = 2
A <- create_A_matrix(1:n, jump_vec, n_seasonal = nbr_sin)
head(A)
A <- create_A_matrix(1:n, jumps = NULL, n_seasonal = nbr_sin)
head(A)

Create a gnssts object

Description

Create a gnssts object

Usage

create.gnssts(t, y, jumps = NULL, sampling_period = 1)
create.gnssts(t, y, jumps = NULL, sampling_period = 1)

Arguments

`t`	A `vector` specifying the time of each observation of the time series.
`y`	A `vector` specifying the values of each observation of the time series.
`jumps`	A `vector` specifying the time values for which there is a jump.
`sampling_period`	An `integer` specifying the sampling period.

Value

A gnssts object.

Examples

phase <- 0.45
amplitude <- 2.5
sigma2_wn <- 15
bias <- 0
trend <- 5 / 365.25
cosU <- amplitude * cos(phase)
sinU <- amplitude * sin(phase)
year <- 5
n <- year * 365
jump_vec <- c(200, 300, 500)
jump_height <- c(10, 15, 20)
nbr_sin <- 1 
A <- create_A_matrix(1:n, jump_vec, n_seasonal = nbr_sin)
x_0 <- c(bias, trend, cosU, sinU, jump_height)
eps <- rnorm(n = n, sd = sqrt(sigma2_wn))
yy <- A %*% x_0 + eps
gnssts_obj <- create.gnssts(t = 1:length(yy), y = yy, jumps = jump_vec)
str(gnssts_obj)

phase <- 0.45
amplitude <- 2.5
sigma2_wn <- 15
bias <- 0
trend <- 5 / 365.25
cosU <- amplitude * cos(phase)
sinU <- amplitude * sin(phase)
year <- 5
n <- year * 365
jump_vec <- c(200, 300, 500)
jump_height <- c(10, 15, 20)
nbr_sin <- 1 
A <- create_A_matrix(1:n, jump_vec, n_seasonal = nbr_sin)
x_0 <- c(bias, trend, cosU, sinU, jump_height)
eps <- rnorm(n = n, sd = sqrt(sigma2_wn))
yy <- A %*% x_0 + eps
gnssts_obj <- create.gnssts(t = 1:length(yy), y = yy, jumps = jump_vec)
str(gnssts_obj)

Estimate a stochastic model in a two-steps procedure using the GMWMX estimator.

Description

Estimate a stochastic model in a two-steps procedure using the GMWMX estimator.

Usage

estimate_gmwmx(
  x,
  theta_0,
  n_seasonal = 1,
  model_string,
  method = "L-BFGS-B",
  maxit = 1e+06,
  ci = FALSE,
  k_iter = 1
)
estimate_gmwmx(
  x,
  theta_0,
  n_seasonal = 1,
  model_string,
  method = "L-BFGS-B",
  maxit = 1e+06,
  ci = FALSE,
  k_iter = 1
)

Arguments

`x`	A `gnssts` object
`theta_0`	A `vector` specifying the initial values for the vector of parameter of the stochastic model considered.
`n_seasonal`	An `integer` specifying the number of seasonal component in the time series.
`model_string`	A `string` specifying the model to be estimated.
`method`	A `string` specifying the numerical optimization method that should be supplied to `optim()`
`maxit`	An `integer` specifying the maximum number of iterations for the numerical optimization procedure.
`ci`	A `boolean` specifying if confidence intervals for the estimated parameters should be computed.
`k_iter`	An `integer` specifying the number of time the two steps GMWMX procedure should be run.

Value

A gnsstsmodel object.

Examples

## Not run: 
data(cola)
fit_gmwmx = estimate_gmwmx(x = cola,
                           theta_0 = c(0.1,0.1,0.1,0.1), 
                           n_seasonal = 1, 
                           model_string = "wn+matern")

## End(Not run)

## Not run: 
data(cola)
fit_gmwmx = estimate_gmwmx(x = cola,
                           theta_0 = c(0.1,0.1,0.1,0.1), 
                           n_seasonal = 1, 
                           model_string = "wn+matern")

## End(Not run)

Estimate a stochastic model based on the MLE and the Hector implementation.

Description

Estimate a stochastic model based on the MLE and the Hector implementation.

Usage

estimate_hector(
  x,
  n_seasonal = 1,
  model_string,
  likelihood_method = "AmmarGrag",
  cleanup = TRUE
)
estimate_hector(
  x,
  n_seasonal = 1,
  model_string,
  likelihood_method = "AmmarGrag",
  cleanup = TRUE
)

Arguments

`x`	A `gnssts` object
`n_seasonal`	An `integer` specifying the number of seasonal component in the time series.
`model_string`	A `string` specifying the model to be estimated.
`likelihood_method`	A `string` taking either value "FullCov" or "AmmarGrag" that specify the method for the Likelihood computation.
`cleanup`	A `boolean` specifying if the files created by the estimation procedure should be cleaned.

Value

A gnsstsmodel object.

Examples

## Not run: 
cola = PBO_get_station(station_name = "COLA", column = "dE", time_range = c(51130, 52000))
fit_mle = estimate_hector(x = cola,
                          n_seasonal = 1, 
                          model_string = "wn+matern")


## End(Not run)

## Not run: 
cola = PBO_get_station(station_name = "COLA", column = "dE", time_range = c(51130, 52000))
fit_mle = estimate_hector(x = cola,
                          n_seasonal = 1, 
                          model_string = "wn+matern")


## End(Not run)

Extract offsets for a PBO station

Description

Extract offsets for a PBO station

Usage

PBO_get_offsets(station_name)
PBO_get_offsets(station_name)

Arguments

station_name

A string specifying the PBO station name.

Value

A vector specifying the offsets of a PBO station.

Examples

## Not run: 
pbo_cola_offsets = PBO_get_offsets(station_name = "COLA")
pbo_cola_offsets

## End(Not run)
## Not run: 
pbo_cola_offsets = PBO_get_offsets(station_name = "COLA")
pbo_cola_offsets

## End(Not run)

Load station data from PBO

Description

Load station data from PBO

Usage

PBO_get_station(station_name, column, time_range = c(-Inf, Inf), scale = 1)
PBO_get_station(station_name, column, time_range = c(-Inf, Inf), scale = 1)

Arguments

`station_name`	A `string` specifying the PBO station name.
`column`	A `string` specifying the name of the column to extract.
`time_range`	A `vector` of 2 specifying the time range of data to extract.
`scale`	A `scalar` specifying an optional scaling parameter applied to the extracted data.

Value

A gnssts object that contains the data associated with the specified PBO station.

Examples

## Not run: 
pbo_cola_data = PBO_get_station("COLA", column="dE")
str(pbo_cola_data)

## End(Not run)
## Not run: 
pbo_cola_data = PBO_get_station("COLA", column="dE")
str(pbo_cola_data)

## End(Not run)

Plotting method for a `gnsstsmodel` object.

Description

Plotting method for a gnsstsmodel object.

Usage

## S3 method for class 'gnsstsmodel'
plot(
  x,
  main = NULL,
  y_unit = "mm",
  x_unit = "days",
  legend_position = "bottomright",
  legend_position_wv = "bottomleft",
  ...
)
## S3 method for class 'gnsstsmodel'
plot(
  x,
  main = NULL,
  y_unit = "mm",
  x_unit = "days",
  legend_position = "bottomright",
  legend_position_wv = "bottomleft",
  ...
)

Arguments

`x`	A `gnsstsmodel` object.
`main`	A `string` specifying the title of the plot.
`y_unit`	A `string` specifying the y axis label of the plot.
`x_unit`	A `string` specifying the x axis label of the plot.
`legend_position`	A `string` specifying the legend position of the plot.
`legend_position_wv`	A `string` specifying the legend position for the wv plot.
`...`	Additional graphical parameters.

Value

No return value. Plot a gnsstsmodel object.

Examples

## Not run: 
data(cola)
fit_gmwmx = estimate_gmwmx(x = cola,
                           theta_0 = c(0.1,0.1,0.1,0.1), 
                           n_seasonal = 1, 
                           model_string = "wn+matern")
plot(fit_gmwmx)

## End(Not run)
## Not run: 
data(cola)
fit_gmwmx = estimate_gmwmx(x = cola,
                           theta_0 = c(0.1,0.1,0.1,0.1), 
                           n_seasonal = 1, 
                           model_string = "wn+matern")
plot(fit_gmwmx)

## End(Not run)

Print method for a `gnsstsmodel` object.

Description

Print method for a gnsstsmodel object.

Usage

## S3 method for class 'gnsstsmodel'
print(x, ...)
## S3 method for class 'gnsstsmodel'
print(x, ...)

Arguments

`x`	A `gnsstsmodel` object.
`...`	Additional graphical parameters.

Value

No return value. Print a gnsstsmodel object.

Examples

## Not run: 
data(cola)
fit_gmwmx = estimate_gmwmx(x = cola,
                           theta_0 = c(0.1,0.1,0.1,0.1), 
                           n_seasonal = 1, 
                           model_string = "wn+matern")
                           
print(fit_gmwmx)

## End(Not run)

## Not run: 
data(cola)
fit_gmwmx = estimate_gmwmx(x = cola,
                           theta_0 = c(0.1,0.1,0.1,0.1), 
                           n_seasonal = 1, 
                           model_string = "wn+matern")
                           
print(fit_gmwmx)

## End(Not run)

Read a gnssts object

Description

Read a gnssts object

Usage

read.gnssts(filename, format = "mom")
read.gnssts(filename, format = "mom")

Arguments

`filename`	A `string` specifying the name of the file to read.
`format`	A `string` specifying the format of the file to read.

Value

Return a gnssts object.

Examples

phase <- 0.45
amplitude <- 2.5
sigma2_wn <- 15
bias <- 0
trend <- 5 / 365.25
cosU <- amplitude * cos(phase)
sinU <- amplitude * sin(phase)
year <- 5
n <- year * 365
jump_vec <- c(200, 300, 500)
jump_height <- c(10, 15, 20)
nbr_sin <- 1
A <- create_A_matrix(1:n, jump_vec, n_seasonal = nbr_sin)
x_0 <- c(bias, trend, cosU, sinU, jump_height)
eps <- rnorm(n = n, sd = sqrt(sigma2_wn))
yy <- A %*% x_0 + eps
gnssts_obj <- create.gnssts(t = 1:length(yy), y = yy, jumps = jump_vec)
str(gnssts_obj)
## Not run: 
write.gnssts(x = gnssts_obj, filename = "test.mom")
gnssts_obj <-read.gnssts(filename = "test.mom", format = "mom")

## End(Not run)
phase <- 0.45
amplitude <- 2.5
sigma2_wn <- 15
bias <- 0
trend <- 5 / 365.25
cosU <- amplitude * cos(phase)
sinU <- amplitude * sin(phase)
year <- 5
n <- year * 365
jump_vec <- c(200, 300, 500)
jump_height <- c(10, 15, 20)
nbr_sin <- 1
A <- create_A_matrix(1:n, jump_vec, n_seasonal = nbr_sin)
x_0 <- c(bias, trend, cosU, sinU, jump_height)
eps <- rnorm(n = n, sd = sqrt(sigma2_wn))
yy <- A %*% x_0 + eps
gnssts_obj <- create.gnssts(t = 1:length(yy), y = yy, jumps = jump_vec)
str(gnssts_obj)
## Not run: 
write.gnssts(x = gnssts_obj, filename = "test.mom")
gnssts_obj <-read.gnssts(filename = "test.mom", format = "mom")

## End(Not run)

Remove outliers from a gnssts object using Hector

Description

Remove outliers from a gnssts object using Hector

Usage

remove_outliers_hector(x, n_seasonal, IQ_factor = 3, cleanup = TRUE)
remove_outliers_hector(x, n_seasonal, IQ_factor = 3, cleanup = TRUE)

Arguments

`x`	A `gnssts` object
`n_seasonal`	An `integer` specifying the number of seasonal component in the time series.
`IQ_factor`	A `double` specifying the number used to scale the interquartile range and corresponding to the argument `IQ_factor` in Hector removeoutliers.ctl
`cleanup`	An `boolean` specifying if temporary files should be cleaned.

Value

A gnssts object.

Examples

phase =     0.45
amplitude = 2.5
sigma2_wn =       15
bias =            0
trend =           5/365.25
cosU =            amplitude*cos(phase)
sinU =            amplitude*sin(phase)
n= 2*365
# define time at which there are jumps
jump_vec =  c(100, 200)
jump_height = c(10, 20)
# generate residuals
eps = rnorm(n = n, sd = sqrt(sigma2_wn))
# add trend, gaps and sin
A = create_A_matrix(1:length(eps), jump_vec, n_seasonal =  1)
# define beta
x_0 = c(bias, trend, jump_height,  cosU,  sinU)
# create time series
yy = A %*% x_0 + eps
plot(yy, type="l")
n_outliers = 30
set.seed(123)
id_outliers=sample(150:350, size = n_outliers)
val_outliers = rnorm(n = n_outliers, mean = max(yy)+10, sd = 5)
yy[id_outliers] = val_outliers
plot(yy, type="l")
# save signal in temp
gnssts_obj = create.gnssts(t = 1:length(yy), y = yy, jumps = jump_vec)
## Not run: 
clean_yy = remove_outliers_hector(x=gnssts_obj, n_seasonal = 1)
plot(clean_yy$t, clean_yy$y, type="l")

## End(Not run)
phase =     0.45
amplitude = 2.5
sigma2_wn =       15
bias =            0
trend =           5/365.25
cosU =            amplitude*cos(phase)
sinU =            amplitude*sin(phase)
n= 2*365
# define time at which there are jumps
jump_vec =  c(100, 200)
jump_height = c(10, 20)
# generate residuals
eps = rnorm(n = n, sd = sqrt(sigma2_wn))
# add trend, gaps and sin
A = create_A_matrix(1:length(eps), jump_vec, n_seasonal =  1)
# define beta
x_0 = c(bias, trend, jump_height,  cosU,  sinU)
# create time series
yy = A %*% x_0 + eps
plot(yy, type="l")
n_outliers = 30
set.seed(123)
id_outliers=sample(150:350, size = n_outliers)
val_outliers = rnorm(n = n_outliers, mean = max(yy)+10, sd = 5)
yy[id_outliers] = val_outliers
plot(yy, type="l")
# save signal in temp
gnssts_obj = create.gnssts(t = 1:length(yy), y = yy, jumps = jump_vec)
## Not run: 
clean_yy = remove_outliers_hector(x=gnssts_obj, n_seasonal = 1)
plot(clean_yy$t, clean_yy$y, type="l")

## End(Not run)

Write a gnssts object

Description

Write a gnssts object

Usage

write.gnssts(x, filename, format = "mom")
write.gnssts(x, filename, format = "mom")

Arguments

`x`	A `R` object to save as a `gnssts` object.
`filename`	A `string` specifying the name of the file to write.
`format`	A `string` specifying the format of the file to write.

Value

No return value. Write a gnssts object in a .mom file by default.

Examples

phase <- 0.45
amplitude <- 2.5
sigma2_wn <- 15
bias <- 0
trend <- 5 / 365.25
cosU <- amplitude * cos(phase)
sinU <- amplitude * sin(phase)
year <- 5
n <- year * 365
jump_vec <- c(200, 300, 500)
jump_height <- c(10, 15, 20)
nbr_sin <- 1
A <- create_A_matrix(1:n, jump_vec, n_seasonal = nbr_sin)
x_0 <- c(bias, trend, cosU, sinU, jump_height)
eps <- rnorm(n = n, sd = sqrt(sigma2_wn))
yy <- A %*% x_0 + eps
gnssts_obj <- create.gnssts(t = 1:length(yy), y = yy, jumps = jump_vec)
str(gnssts_obj)
## Not run: 
write.gnssts(x = gnssts_obj, filename = "test.mom")

## End(Not run)
phase <- 0.45
amplitude <- 2.5
sigma2_wn <- 15
bias <- 0
trend <- 5 / 365.25
cosU <- amplitude * cos(phase)
sinU <- amplitude * sin(phase)
year <- 5
n <- year * 365
jump_vec <- c(200, 300, 500)
jump_height <- c(10, 15, 20)
nbr_sin <- 1
A <- create_A_matrix(1:n, jump_vec, n_seasonal = nbr_sin)
x_0 <- c(bias, trend, cosU, sinU, jump_height)
eps <- rnorm(n = n, sd = sqrt(sigma2_wn))
yy <- A %*% x_0 + eps
gnssts_obj <- create.gnssts(t = 1:length(yy), y = yy, jumps = jump_vec)
str(gnssts_obj)
## Not run: 
write.gnssts(x = gnssts_obj, filename = "test.mom")

## End(Not run)

Package 'gmwmx'

Help Index

GNSS time series for PBO station COLA

Description

Usage

Format

Source

Compare graphically two gnsstsmodel objects.

Description

Usage

Arguments

Value

Examples

Define matrix A of the functional model

Description

Usage

Arguments

Value

Examples

Create a gnssts object

Description

Usage

Arguments

Value

Examples

Estimate a stochastic model in a two-steps procedure using the GMWMX estimator.

Description

Usage

Arguments

Value

Examples

Estimate a stochastic model based on the MLE and the Hector implementation.

Description

Usage

Arguments

Value

Examples

Extract offsets for a PBO station

Description

Usage

Arguments

Value

Examples

Load station data from PBO

Description

Usage

Arguments

Value

Examples

Plotting method for a gnsstsmodel object.

Description

Usage

Arguments

Value

Examples

Print method for a gnsstsmodel object.

Description

Usage

Arguments

Value

Examples

Read a gnssts object

Description

Usage

Arguments

Value

Examples

Remove outliers from a gnssts object using Hector

Description

Usage

Arguments

Value

Examples

Write a gnssts object

Description

Usage

Arguments

Value

Examples

Compare graphically two `gnsstsmodel` objects.

Plotting method for a `gnsstsmodel` object.

Print method for a `gnsstsmodel` object.