Title: | Time Series Analysis with State Space Model |
---|---|
Description: | Functions for statistical analysis, modeling and simulation of time series with state space model, based on the methodology in Kitagawa (2020, ISBN: 978-0-367-18733-0). |
Authors: | The Institute of Statistical Mathematics, based on the program by Genshiro Kitagawa |
Maintainer: | Masami Saga <[email protected]> |
License: | GPL (>= 2) |
Version: | 1.3.4-5 |
Built: | 2024-12-10 06:36:58 UTC |
Source: | CRAN |
R functions for statistical analysis, modeling and simulation of time series with state space model.
This package provides functions for statistical analysis, modeling and simulation of time series. These functions are developed based on source code of "FORTRAN 77 Programming for Time Series Analysis".
After that, the revised edition "Introduction to Time Series Analysis (in Japanese)" and the translation version "Introduction to Time Series Modeling" are published.
Currently the revised edition "Introduction to Time Series Modeling with Applications in R" is published, in which calculations of most of the modeling or methods are explained using this package.
Kitagawa, G. and Gersch, W. (1996) Smoothness Priors Analysis of Time Series. Lecture Notes in Statistics, No.116, Springer-Verlag.
Kitagawa, G. (2010) Introduction to Time Series Modeling. Chapman & Hall/CRC.
Kitagawa, G. (2020) Introduction to Time Series Modeling with Applications in R. Chapman & Hall/CRC.
Kitagawa, G. (1993) FORTRAN 77 Programming for Time Series Analysis. The Iwanami Computer Science Series, Iwanami Publishing Company (in Japanese).
Kitagawa, G. (2005) Introduction to Time Series Analysis. Iwanami Publishing Company (in Japanese).
Kitagawa, G. (2020) Introduction to Time Series Modeling with R. Iwanami Publishing Company (in Japanese).
Fit a univariate AR model by the Yule-Walker method, the least squares (Householder) method or the PARCOR method.
arfit(y, lag = NULL, method = 1, plot = TRUE, ...)
arfit(y, lag = NULL, method = 1, plot = TRUE, ...)
y |
a univariate time series. |
||||||||||
lag |
highest order of AR model. Default is |
||||||||||
method |
estimation procedure.
|
||||||||||
plot |
logical. If |
||||||||||
... |
graphical arguments passed to the |
An object of class "arfit"
which has a plot
method. This is a
list with the following components:
sigma2 |
innovation variance. |
maice.order |
order of minimum AIC. |
aic |
AICs of the estimated AR models. |
arcoef |
AR coefficients of the estimated AR models. |
parcor |
PARCOR. |
spec |
power spectrum (in log scale) of the AIC best AR model. |
tsname |
the name of the univariate time series |
Kitagawa, G. (2020) Introduction to Time Series Modeling with Applications in R. Chapman & Hall/CRC.
# Sunspot number data data(Sunspot) arfit(log10(Sunspot), lag = 20, method = 1) # BLSALLFOOD data data(BLSALLFOOD) arfit(BLSALLFOOD)
# Sunspot number data data(Sunspot) arfit(log10(Sunspot), lag = 20, method = 1) # BLSALLFOOD data data(BLSALLFOOD) arfit(BLSALLFOOD)
Calculate impulse response function, autocovariance function, autocorrelation function and characteristic roots of given scalar ARMA model.
armachar(arcoef = NULL, macoef = NULL, v, lag = 50, nf = 200, plot = TRUE, ...)
armachar(arcoef = NULL, macoef = NULL, v, lag = 50, nf = 200, plot = TRUE, ...)
arcoef |
AR coefficients. |
macoef |
MA coefficients. |
v |
innovation variance. |
lag |
maximum lag of autocovariance function. |
nf |
number of frequencies in evaluating spectrum. |
plot |
logical. If |
... |
graphical arguments passed to the |
The ARMA model is given by
where is AR order,
is MA order and
is a zero
mean white noise.
Characteristic roots of AR / MA operator is a list with the following components:
re: real part
im: imaginary part
amp:
atan:
degree
An object of class "arma"
which has a plot
method. This is a
list with components:
impuls |
impulse response function. |
acov |
autocovariance function. |
parcor |
PARCOR. |
spec |
power spectrum. |
croot.ar |
characteristic roots of AR operator. See Details. |
croot.ma |
characteristic roots of MA operator. See Details. |
Kitagawa, G. (2020) Introduction to Time Series Modeling with Applications in R. Chapman & Hall/CRC.
# AR model : y(n) = a(1)*y(n-1) + a(2)*y(n-2) + v(n) a <- c(0.9 * sqrt(3), -0.81) armachar(arcoef = a, v = 1.0, lag = 20) # MA model : y(n) = v(n) - b(1)*v(n-1) - b(2)*v(n-2) b <- c(0.9 * sqrt(2), -0.81) armachar(macoef = b, v = 1.0, lag = 20) # ARMA model : y(n) = a(1)*y(n-1) + a(2)*y(n-2) # + v(n) - b(1)*v(n-1) - b(2)*v(n-2) armachar(arcoef = a, macoef = b, v = 1.0, lag = 20)
# AR model : y(n) = a(1)*y(n-1) + a(2)*y(n-2) + v(n) a <- c(0.9 * sqrt(3), -0.81) armachar(arcoef = a, v = 1.0, lag = 20) # MA model : y(n) = v(n) - b(1)*v(n-1) - b(2)*v(n-2) b <- c(0.9 * sqrt(2), -0.81) armachar(macoef = b, v = 1.0, lag = 20) # ARMA model : y(n) = a(1)*y(n-1) + a(2)*y(n-2) # + v(n) - b(1)*v(n-1) - b(2)*v(n-2) armachar(arcoef = a, macoef = b, v = 1.0, lag = 20)
Fit a scalar ARMA model by maximum likelihood method.
armafit(y, ar.order, ar = NULL, ma.order, ma = NULL)
armafit(y, ar.order, ar = NULL, ma.order, ma = NULL)
y |
a univariate time series. |
ar.order |
AR order. |
ar |
initial AR coefficients. If |
ma.order |
MA order. |
ma |
initial MA coefficients. If |
sigma2 |
innovation variance. |
llkhood |
log-likelihood of the model. |
aic |
AIC of the model. |
arcoef |
AR coefficients. |
macoef |
MA coefficients. |
Kitagawa, G. (2020) Introduction to Time Series Modeling with Applications in R. Chapman & Hall/CRC.
# Sunspot number data data(Sunspot) y <- log10(Sunspot) z <- armafit(y, ar.order = 3, ma.order = 3) z armachar(arcoef = z$arcoef, macoef = z$macoef, v = z$sigma2, lag = 20)
# Sunspot number data data(Sunspot) y <- log10(Sunspot) z <- armafit(y, ar.order = 3, ma.order = 3) z armachar(arcoef = z$arcoef, macoef = z$macoef, v = z$sigma2, lag = 20)
Estimate all ARMA models within the user-specified maximum order by maximum likelihood method.
armafit2(y, ar.order, ma.order)
armafit2(y, ar.order, ma.order)
y |
a univariate time series. |
ar.order |
maximum AR order. |
ma.order |
maximum MA order. |
aicmin |
minimum AIC. |
maice.order |
AR and MA orders of minimum AIC model. |
sigma2 |
innovation variance of all models. |
llkhood |
log-likelihood of all models. |
aic |
AIC of all models. |
coef |
AR and MA coefficients of all models. |
Kitagawa, G. (2020) Introduction to Time Series Modeling with Applications in R. Chapman & Hall/CRC.
# Sunspot number data data(Sunspot) y <- log10(Sunspot) armafit2(y, ar.order = 5, ma.order = 5)
# Sunspot number data data(Sunspot) y <- log10(Sunspot) armafit2(y, ar.order = 5, ma.order = 5)
The monthly time series of the number of workers engaged in food industries in the United States (January 1967 - December 1979).
data(BLSALLFOOD)
data(BLSALLFOOD)
A time series of 156 observations.
The data were obtained from the United States Bureau of Labor Statistics (BLS).
Compute Box-Cox transformation and find an optimal lambda with minimum AIC.
boxcox(y, plot = TRUE, ...)
boxcox(y, plot = TRUE, ...)
y |
a univariate time series. |
plot |
logical. If |
... |
graphical arguments passed to |
An object of class "boxcox"
, which is a list with the following
components:
mean |
mean of original data. |
var |
variance of original data. |
aic |
AIC of the model with respect to the original data. |
llkhood |
log-likelihood of the model with respect to the original data. |
z |
transformed data with the AIC best lambda. |
aic.z |
AIC of the model with respect to the transformed data. |
llkhood.z |
log-likelihood of the model with respect to the transformed data. |
Kitagawa, G. (2020) Introduction to Time Series Modeling with Applications in R. Chapman & Hall/CRC.
# Sunspot number data data(Sunspot) boxcox(Sunspot) # Wholesale hardware data data(WHARD) boxcox(WHARD)
# Sunspot number data data(Sunspot) boxcox(Sunspot) # Wholesale hardware data data(WHARD) boxcox(WHARD)
Compute cross-covariance and cross-correlation functions of the multivariate time series.
crscor(y, lag = NULL, outmin = NULL, outmax = NULL, plot = TRUE, ...)
crscor(y, lag = NULL, outmin = NULL, outmax = NULL, plot = TRUE, ...)
y |
a multivariate time series. |
lag |
maximum lag. Default is |
outmin |
bound for outliers in low side. A default value is -1.0e+30 for each dimension. |
outmax |
bound for outliers in high side. A default value is 1.0e+30 for each dimension. |
plot |
logical. If |
... |
graphical arguments passed to the |
An object of class "crscor"
which has a plot
method. This is a
list with the following components:
cov |
cross-covariances. |
cor |
cross-correlations. |
mean |
mean vector. |
Kitagawa, G. (2020) Introduction to Time Series Modeling with Applications in R. Chapman & Hall/CRC.
# Yaw rate, rolling, pitching and rudder angle of a ship data(HAKUSAN) y <- as.matrix(HAKUSAN[, 2:4]) # Rolling, Pitching, Rudder crscor(y, lag = 50) # The groundwater level and the atmospheric pressure data(Haibara) crscor(Haibara, lag = 50)
# Yaw rate, rolling, pitching and rudder angle of a ship data(HAKUSAN) y <- as.matrix(HAKUSAN[, 2:4]) # Rolling, Pitching, Rudder crscor(y, lag = 50) # The groundwater level and the atmospheric pressure data(Haibara) crscor(Haibara, lag = 50)
Compute a periodogram of the univariate time series via FFT.
fftper(y, window = 1, plot = TRUE, ...)
fftper(y, window = 1, plot = TRUE, ...)
y |
a univariate time series. |
window |
smoothing window type. (0: box-car, 1: Hanning, 2: Hamming) |
plot |
logical. If |
... |
graphical arguments passed to |
Hanning Window : | = 0.5 |
= 0.25 |
Hamming Window : | = 0.54 |
= 0.23 |
An object of class "spg"
, which is a list with the following components:
period |
periodogram. |
smoothed.period |
smoothed periodogram. If there is not a negative number, logarithm of smoothed periodogram. |
log.scale |
logical. If |
tsname |
the name of the univariate time series |
We assume that the length of the input time series
y
is a power
of 2. If is not a power of 2, calculate using the FFT by appending 0's
behind the data
y
.
Kitagawa, G. (2020) Introduction to Time Series Modeling with Applications in R. Chapman & Hall/CRC.
# Yaw rate, rolling, pitching and rudder angle of a ship data(HAKUSAN) YawRate <- HAKUSAN[, 1] fftper(YawRate, window = 0)
# Yaw rate, rolling, pitching and rudder angle of a ship data(HAKUSAN) YawRate <- HAKUSAN[, 1] fftper(YawRate, window = 0)
A bivariate time series of the groundwater level and the atmospheric pressure that were observed at 10-minuite intervals at the Haibara observatory of the Tokai region, Japan.
data(Haibara)
data(Haibara)
A data frame with 400 observations on the following 2 variables.
[, 1] | Groundwater level |
[, 2] | Atmospheric pressure |
The data were offered by Dr. M. Takahashi and Dr. N. Matsumoto of National Institute of Advanced Industrial Science and Technology.
data(Haibara) ## put histograms on the diagonal panel.hist <- function(x, ...) { usr <- par("usr") par(usr = c(usr[1:2], 0, 1.3)) nB <- 15; nB1 <- nB + 1 xmin <- min(x, na.rm = TRUE) xmax <- max(x, na.rm = TRUE) w <- (xmax - xmin) / nB breaks <- xmin b <- xmin for (i in 1:nB) { b <- b + w breaks <- c(breaks, b) } h <- hist(x, breaks = breaks, plot = FALSE) y <- h$counts y <- y / max(y) rect(breaks[1:nB], 0, breaks[2:nB1], y, ...) } par(xaxs = "i", yaxs = "i", xaxt = "n", yaxt = "n") pairs(Haibara, diag.panel = panel.hist, pch = 20, cex.labels = 1.5, label.pos = 0.9, lower.panel = NULL)
data(Haibara) ## put histograms on the diagonal panel.hist <- function(x, ...) { usr <- par("usr") par(usr = c(usr[1:2], 0, 1.3)) nB <- 15; nB1 <- nB + 1 xmin <- min(x, na.rm = TRUE) xmax <- max(x, na.rm = TRUE) w <- (xmax - xmin) / nB breaks <- xmin b <- xmin for (i in 1:nB) { b <- b + w breaks <- c(breaks, b) } h <- hist(x, breaks = breaks, plot = FALSE) y <- h$counts y <- y / max(y) rect(breaks[1:nB], 0, breaks[2:nB1], y, ...) } par(xaxs = "i", yaxs = "i", xaxt = "n", yaxt = "n") pairs(Haibara, diag.panel = panel.hist, pch = 20, cex.labels = 1.5, label.pos = 0.9, lower.panel = NULL)
A multivariate time series of a ship's yaw rate, rolling, pitching and rudder angles which were recorded every second while navigating across the Pacific Ocean.
data(HAKUSAN)
data(HAKUSAN)
A data frame with 1000 observations on the following 4 variables.
[, 1] | YawRate | yaw rate |
[, 2] | Rolling | rolling |
[, 3] | Pitching | pitching |
[, 4] | Rudder | rudder angle |
The data were offered by Prof. K. Ohtsu of Tokyo University of Marine Science and Technology.
data(HAKUSAN) HAKUSAN234 <- HAKUSAN[, c(2,3,4)] ## put histograms on the diagonal panel.hist <- function(x, ...) { usr <- par("usr") par(usr = c(usr[1:2], 0, 1.3)) nB <- 20; nB1 <- nB + 1 xmin <- min(x) xmax <- max(x) w <- (xmax - xmin) / nB breaks <- xmin b <- xmin for (i in 1:nB) { b <- b + w breaks <- c(breaks, b) } h <- hist(x, breaks = breaks, plot = FALSE) y <- h$counts; y <- y / max(y) rect(breaks[1:nB], 0, breaks[2:nB1], y, ...) } par(xaxs = "i", yaxs = "i", xaxt = "n", yaxt = "n") pairs(HAKUSAN234, diag.panel = panel.hist, pch = 20, cex.labels = 1.5, label.pos = 0.9, lower.panel = NULL)
data(HAKUSAN) HAKUSAN234 <- HAKUSAN[, c(2,3,4)] ## put histograms on the diagonal panel.hist <- function(x, ...) { usr <- par("usr") par(usr = c(usr[1:2], 0, 1.3)) nB <- 20; nB1 <- nB + 1 xmin <- min(x) xmax <- max(x) w <- (xmax - xmin) / nB breaks <- xmin b <- xmin for (i in 1:nB) { b <- b + w breaks <- c(breaks, b) } h <- hist(x, breaks = breaks, plot = FALSE) y <- h$counts; y <- y / max(y) rect(breaks[1:nB], 0, breaks[2:nB1], y, ...) } par(xaxs = "i", yaxs = "i", xaxt = "n", yaxt = "n") pairs(HAKUSAN234, diag.panel = panel.hist, pch = 20, cex.labels = 1.5, label.pos = 0.9, lower.panel = NULL)
Compute Kullback-Leibler information.
klinfo(distg = 1, paramg = c(0, 1), distf = 1, paramf, xmax = 10)
klinfo(distg = 1, paramg = c(0, 1), distf = 1, paramf, xmax = 10)
distg |
function for the true density (1 or 2).
|
||||||||||||
paramg |
parameter vector of true density. |
||||||||||||
distf |
function for the model density (1 or 2).
|
||||||||||||
paramf |
parameter vector of the model density. |
||||||||||||
xmax |
upper limit of integration. lower limit xmin = - |
nint |
number of function evaluation. |
dx |
delta. |
KLI |
Kullback-Leibler information, |
gint |
integration of |
Kitagawa, G. (2020) Introduction to Time Series Modeling with Applications in R. Chapman & Hall/CRC.
# g:Gauss, f:Gauss klinfo(distg = 1, paramg = c(0, 1), distf = 1, paramf = c(0.1, 1.5), xmax = 8) # g:Gauss, f:Cauchy klinfo(distg = 1, paramg = c(0, 1), distf = 2, paramf = c(0, 1), xmax = 8)
# g:Gauss, f:Gauss klinfo(distg = 1, paramg = c(0, 1), distf = 1, paramf = c(0.1, 1.5), xmax = 8) # g:Gauss, f:Cauchy klinfo(distg = 1, paramg = c(0, 1), distf = 2, paramf = c(0, 1), xmax = 8)
Decompose time series to stationary subintervals and estimate local spectrum.
lsar(y, max.arorder = 20, ns0, plot = TRUE, ...)
lsar(y, max.arorder = 20, ns0, plot = TRUE, ...)
y |
a univariate time series. |
max.arorder |
highest order of AR model. |
ns0 |
length of basic local span. |
plot |
logical. If |
... |
graphical arguments passed to the |
An object of class "lsar"
which has a plot
method. This is a
list with the following components:
model |
1: pooled model is accepted. |
ns |
number of observations of local span. |
span |
start points and end points of local spans. |
nf |
number of frequencies in computing local power spectrum. |
ms |
order of switched model. |
sds |
innovation variance of switched model. |
aics |
AIC of switched model. |
mp |
order of pooled model. |
sdp |
innovation variance of pooled model. |
aics |
AIC of pooled model. |
spec |
local spectrum. |
tsname |
the name of the univariate time series |
Kitagawa, G. (2020) Introduction to Time Series Modeling with Applications in R. Chapman & Hall/CRC.
# seismic data data(MYE1F) lsar(MYE1F, max.arorder = 10, ns0 = 100)
# seismic data data(MYE1F) lsar(MYE1F, max.arorder = 10, ns0 = 100)
Precisely estimate a change point of subinterval for locally stationary AR model.
lsar.chgpt(y, max.arorder = 20, subinterval, candidate, plot = TRUE, ...)
lsar.chgpt(y, max.arorder = 20, subinterval, candidate, plot = TRUE, ...)
y |
a univariate time series. |
max.arorder |
highest order of AR model. |
subinterval |
a vector of the form |
candidate |
a vector of the form
|
plot |
logical. If |
... |
graphical arguments passed to the |
An object of class "chgpt"
which has a plot
method. This is a
list with the following components:
aic |
AICs of the AR models fitted on |
aicmin |
minimum AIC. |
change.point |
estimated change point. |
subint |
information about the original sub-interval. |
Kitagawa, G. (2020) Introduction to Time Series Modeling with Applications in R. Chapman & Hall/CRC.
# seismic data data(MYE1F) lsar.chgpt(MYE1F, max.arorder = 10, subinterval = c(200, 1000), candidate = c(400, 800)) lsar.chgpt(MYE1F, max.arorder = 10, subinterval = c(600, 1400), candidate = c(800, 1200))
# seismic data data(MYE1F) lsar.chgpt(MYE1F, max.arorder = 10, subinterval = c(200, 1000), candidate = c(400, 800)) lsar.chgpt(MYE1F, max.arorder = 10, subinterval = c(600, 1400), candidate = c(800, 1200))
Compute regression coefficients of the model with minimum AIC by the least squares method via Householder transformation.
lsqr(y, lag = NULL, period = 365, plot = TRUE, ...)
lsqr(y, lag = NULL, period = 365, plot = TRUE, ...)
y |
a univariate time series. |
lag |
number of sine and cosine components. Default is |
period |
period of one cycle. |
plot |
logical. If |
... |
graphical arguments passed to |
An object of class "lsqr"
, which is a list with the following
components:
aic |
AIC's of the model with order |
sigma2 |
residual variance of the model with order |
maice.order |
order of minimum AIC. |
regress |
regression coefficients of the model. |
tripoly |
trigonometric polynomial. |
Kitagawa, G. (2020) Introduction to Time Series Modeling with Applications in R. Chapman & Hall/CRC.
# The daily maximum temperatures in Tokyo data(Temperature) lsqr(Temperature, lag = 10)
# The daily maximum temperatures in Tokyo data(Temperature) lsqr(Temperature, lag = 10)
Fit a multivariate AR model by the Yule-Walker method.
marfit(y, lag = NULL)
marfit(y, lag = NULL)
y |
a multivariate time series. |
lag |
highest order of fitted AR models. Default is |
An object of class "maryule"
, which is a list with the following
components:
maice.order |
order of minimum AIC. |
aic |
AIC's of the AR models with order |
v |
innovation covariance matrix of the AIC best model. |
arcoef |
AR coefficients of the AIC best model. |
Kitagawa, G. (2020) Introduction to Time Series Modeling with Applications in R. Chapman & Hall/CRC.
# Yaw rate, rolling, pitching and rudder angle of a ship data(HAKUSAN) yy <- as.matrix(HAKUSAN[, c(1,2,4)]) # Yaw rate, Pitching, Rudder angle nc <- dim(yy)[1] n <- seq(1, nc, by = 2) y <- yy[n, ] marfit(y, 20)
# Yaw rate, rolling, pitching and rudder angle of a ship data(HAKUSAN) yy <- as.matrix(HAKUSAN[, c(1,2,4)]) # Yaw rate, Pitching, Rudder angle nc <- dim(yy)[1] n <- seq(1, nc, by = 2) y <- yy[n, ] marfit(y, 20)
Fit a multivariate AR model by least squares method.
marlsq(y, lag = NULL)
marlsq(y, lag = NULL)
y |
a multivariate time series. |
lag |
highest AR order. Default is |
An object of class "marlsq"
, which is a list with the following
components:
maice.order |
order of the MAICE model. |
aic |
AIC of the MAR model with minimum AIC orders. |
v |
innovation covariance matrix. |
arcoef |
AR coefficient matrices. |
Kitagawa, G. (2020) Introduction to Time Series Modeling with Applications in R. Chapman & Hall/CRC.
# Yaw rate, rolling, pitching and rudder angle of a ship data(HAKUSAN) y <- as.matrix(HAKUSAN[, c(1,2,4)]) # Yaw rate, Rolling, Rudder angle z <- marlsq(y) z marspc(z$arcoef, v = z$v)
# Yaw rate, rolling, pitching and rudder angle of a ship data(HAKUSAN) y <- as.matrix(HAKUSAN[, c(1,2,4)]) # Yaw rate, Rolling, Rudder angle z <- marlsq(y) z marspc(z$arcoef, v = z$v)
Compute cross spectra, coherency and power contribution.
marspc(arcoef, v, plot = TRUE, ...)
marspc(arcoef, v, plot = TRUE, ...)
arcoef |
AR coefficient matrices. |
v |
innovation variance matrix. |
plot |
logical. If |
... |
graphical arguments passed to the |
An object of class "marspc"
which has a plot
method. This is a
list with the following components:
spec |
cross spectra. |
amp |
amplitude spectra. |
phase |
phase spectra. |
coh |
simple coherency. |
power |
decomposition of power spectra. |
rpower |
relative power contribution. |
Kitagawa, G. (2020) Introduction to Time Series Modeling with Applications in R. Chapman & Hall/CRC.
# Yaw rate, rolling, pitching and rudder angle of a ship data(HAKUSAN) yy <- as.matrix(HAKUSAN[, c(1,2,4)]) nc <- dim(yy)[1] n <- seq(1, nc, by = 2) y <- yy[n, ] z <- marfit(y, lag = 20) marspc(z$arcoef, v = z$v)
# Yaw rate, rolling, pitching and rudder angle of a ship data(HAKUSAN) yy <- as.matrix(HAKUSAN[, c(1,2,4)]) nc <- dim(yy)[1] n <- seq(1, nc, by = 2) y <- yy[n, ] z <- marfit(y, lag = 20) marspc(z$arcoef, v = z$v)
The time series of East-West components of seismic waves, recorded every 0.02 seconds.
data(MYE1F)
data(MYE1F)
A time series of 2600 observations.
Takanami, T. (1991), "ISM data 43-3-01: Seismograms of foreshocks of 1982 Urakawa-Oki earthquake", Ann. Inst. Statist. Math., 43, 605.
Simulation by non-Gaussian state space model.
ngsim(n = 200, trend = NULL, seasonal.order = 0, seasonal = NULL, arcoef = NULL, ar = NULL, noisew = 1, wminmax = NULL, paramw = NULL, noisev = 1, vminmax = NULL, paramv = NULL, seed = NULL, plot = TRUE, ...)
ngsim(n = 200, trend = NULL, seasonal.order = 0, seasonal = NULL, arcoef = NULL, ar = NULL, noisew = 1, wminmax = NULL, paramw = NULL, noisev = 1, vminmax = NULL, paramv = NULL, seed = NULL, plot = TRUE, ...)
n |
number of data generated by simulation. |
||||||||||||||
trend |
initial values of trend component of length |
||||||||||||||
seasonal.order |
order of seasonal component model (0, 1, 2). |
||||||||||||||
seasonal |
if |
||||||||||||||
arcoef |
AR coefficients. |
||||||||||||||
ar |
initial values of AR component. |
||||||||||||||
noisew |
type of the observational noise.
|
||||||||||||||
wminmax |
lower and upper bound of observational noise. |
||||||||||||||
paramw |
parameter of the observational noise density.
|
||||||||||||||
noisev |
type of the system noise.
|
||||||||||||||
vminmax |
lower and upper bound of system noise. |
||||||||||||||
paramv |
parameter of the system noise density.
|
||||||||||||||
seed |
arbitrary positive integer to generate a sequence of uniform random numbers. The default seed is based on the current time. |
||||||||||||||
plot |
logical. If |
||||||||||||||
... |
graphical arguments passed to |
An object of class "simulate"
, giving simulated data of non-Gaussian
state space model.
Kitagawa, G. (2020) Introduction to Time Series Modeling with Applications in R. Chapman & Hall/CRC.
ar1 <- ngsim(n = 400, arcoef = 0.95, noisew = 1, paramw = 1, noisev = 1, paramv = 1, seed = 555) plot(ar1, use = c(201, 400)) ar2 <- ngsim(n = 400, arcoef = c(1.3, -0.8), noisew = 1, paramw = 1, noisev = 1, paramv = 1, seed = 555) plot(ar2, use = c(201, 400))
ar1 <- ngsim(n = 400, arcoef = 0.95, noisew = 1, paramw = 1, noisev = 1, paramv = 1, seed = 555) plot(ar1, use = c(201, 400)) ar2 <- ngsim(n = 400, arcoef = c(1.3, -0.8), noisew = 1, paramw = 1, noisev = 1, paramv = 1, seed = 555) plot(ar2, use = c(201, 400))
Trend estimation by non-Gaussian smoothing.
ngsmth(y, noisev = 2, tau2, bv = 1.0, noisew = 1, sigma2, bw = 1.0, initd = 1, k = 200, plot = TRUE, ...)
ngsmth(y, noisev = 2, tau2, bv = 1.0, noisew = 1, sigma2, bw = 1.0, initd = 1, k = 200, plot = TRUE, ...)
y |
a univariate time series. |
||||||||
noisev |
type of system noise density.
|
||||||||
tau2 |
variance or dispersion of system noise. |
||||||||
bv |
shape parameter of system noise (for |
||||||||
noisew |
type of observation noise density
|
||||||||
sigma2 |
variance or dispersion of observation noise. |
||||||||
bw |
shape parameter of observation noise (for |
||||||||
initd |
type of density function.
|
||||||||
k |
number of intervals in numerical integration. |
||||||||
plot |
logical. If |
||||||||
... |
graphical arguments passed to |
Consider a one-dimensional state space model
where the observation noise is assumed to be Gaussian
distributed and the system noise
is assumed to be distributed
as the Pearson system
with and
.
This broad family of distributions includes the Cauchy distribution
() and
-distribution (
).
An object of class "ngsmth"
, which is a list with the following
components:
llkhood |
log-likelihood. |
trend |
trend. |
smt |
smoothed density. |
Kitagawa, G. (2020) Introduction to Time Series Modeling with Applications in R. Chapman & Hall/CRC.
Kitagawa, G. and Gersch, W. (1996) Smoothness Priors Analysis of Time Series. Lecture Notes in Statistics, No.116, Springer-Verlag.
## test data data(PfilterSample) par(mar = c(3, 3, 1, 1) + 0.1) # system noise density : Gaussian (normal) s1 <- ngsmth(PfilterSample, noisev = 1, tau2 = 1.4e-02, noisew = 1, sigma2 = 1.048) s1 plot(s1, "smt", theta = 25, phi = 30, expand = 0.25, col = "white") # system noise density : Pearson family s2 <- ngsmth(PfilterSample, noisev = 2, tau2 = 2.11e-10, bv = 0.6, noisew = 1, sigma2 = 1.042) s2 plot(s2, "smt", theta = 25, phi = 30, expand = 0.25, col = "white") ## seismic data data(MYE1F) n <- length(MYE1F) yy <- rep(0, n) for (i in 2:n) yy[i] <- MYE1F[i] - 0.5 * MYE1F[i-1] m <- seq(1, n, by = 2) y <- yy[m] z <- tvvar(y, trend.order = 2, tau2.ini = 4.909e-02, delta = 1.0e-06) # system noise density : Gaussian (normal) s3 <- ngsmth(z$sm, noisev = 1, tau2 = z$tau2, noisew = 2, sigma2 = pi*pi/6, k = 190) s3 plot(s3, "smt", phi = 50, expand = 0.5, col = 8)
## test data data(PfilterSample) par(mar = c(3, 3, 1, 1) + 0.1) # system noise density : Gaussian (normal) s1 <- ngsmth(PfilterSample, noisev = 1, tau2 = 1.4e-02, noisew = 1, sigma2 = 1.048) s1 plot(s1, "smt", theta = 25, phi = 30, expand = 0.25, col = "white") # system noise density : Pearson family s2 <- ngsmth(PfilterSample, noisev = 2, tau2 = 2.11e-10, bv = 0.6, noisew = 1, sigma2 = 1.042) s2 plot(s2, "smt", theta = 25, phi = 30, expand = 0.25, col = "white") ## seismic data data(MYE1F) n <- length(MYE1F) yy <- rep(0, n) for (i in 2:n) yy[i] <- MYE1F[i] - 0.5 * MYE1F[i-1] m <- seq(1, n, by = 2) y <- yy[m] z <- tvvar(y, trend.order = 2, tau2.ini = 4.909e-02, delta = 1.0e-06) # system noise density : Gaussian (normal) s3 <- ngsmth(z$sm, noisev = 1, tau2 = z$tau2, noisew = 2, sigma2 = pi*pi/6, k = 190) s3 plot(s3, "smt", phi = 50, expand = 0.5, col = 8)
A daily closing values of the Japanese stock price index, Nikkei225, quoted from January 4, 1988, to December 30, 1993.
data(Nikkei225)
data(Nikkei225)
A time series of 1480 observations.
https://indexes.nikkei.co.jp/nkave/archives/data
The series generated by the nonlinear state-space model.
data(NLmodel)
data(NLmodel)
A matrix with 100 rows and 2 columns.
[, 1] | |
[, 2] |
|
The system model and the observation model
are
generated by following state-space model:
where ~
,
~
,
~
.
Evaluate probability density function for normal distribution, Cauchy distribution, Pearson distribution, exponential distribution, Chi-square distributions, double exponential distribution and uniform distribution.
pdfunc(model = "norm", mean = 0, sigma2 = 1, mu = 0, tau2 = 1, shape, lambda = 1, side = 1, df, xmin = 0, xmax = 1, plot = TRUE, ...)
pdfunc(model = "norm", mean = 0, sigma2 = 1, mu = 0, tau2 = 1, shape, lambda = 1, side = 1, df, xmin = 0, xmax = 1, plot = TRUE, ...)
model |
a character string indicating the model type of probability
density function: either |
mean |
mean. (valid for |
sigma2 |
variance. (valid for |
mu |
location parameter |
tau2 |
dispersion parameter |
shape |
shape parameter (> 0.5). (valid for |
lambda |
lambda |
side |
1: exponential, 2: two-sided exponential.
(valid for |
df |
degree of freedoms |
xmin |
lower bound of the interval. |
xmax |
upper bound of the interval. |
plot |
logical. If |
... |
graphical arguments passed to the |
An object of class "pdfunc"
which has a plot
method. This is a
list with the following components:
density |
values of density function. |
interval |
lower and upper bound of interval. |
param |
parameters of model. |
Kitagawa, G. (2020) Introduction to Time Series Modeling with Applications in R. Chapman & Hall/CRC.
# normal distribution pdfunc(model = "norm", xmin = -4, xmax = 4) # Cauchy distribution pdfunc(model = "Cauchy", xmin = -4, xmax = 4) # Pearson distribution pdfunc(model = "Pearson", shape = 2, xmin = -4, xmax = 4) # exponential distribution pdfunc(model = "exp", xmin = 0, xmax = 8) pdfunc(model = "exp", xmin = -4, xmax = 4) # Chi-square distribution pdfunc(model = "Chi2", df = 3, xmin = 0, xmax = 8) # double exponential distribution pdfunc(model = "dexp", xmin = -4, xmax = 2) # uniform distribution pdfunc(model = "unif", xmin = 0, xmax = 1)
# normal distribution pdfunc(model = "norm", xmin = -4, xmax = 4) # Cauchy distribution pdfunc(model = "Cauchy", xmin = -4, xmax = 4) # Pearson distribution pdfunc(model = "Pearson", shape = 2, xmin = -4, xmax = 4) # exponential distribution pdfunc(model = "exp", xmin = 0, xmax = 8) pdfunc(model = "exp", xmin = -4, xmax = 4) # Chi-square distribution pdfunc(model = "Chi2", df = 3, xmin = 0, xmax = 8) # double exponential distribution pdfunc(model = "dexp", xmin = -4, xmax = 2) # uniform distribution pdfunc(model = "unif", xmin = 0, xmax = 1)
Compute a periodogram of the univariate time series.
period(y, window = 1, lag = NULL, minmax = c(-1.0e+30, 1.0e+30), plot = TRUE, ...)
period(y, window = 1, lag = NULL, minmax = c(-1.0e+30, 1.0e+30), plot = TRUE, ...)
y |
a univariate time series. |
window |
smoothing window type. (0: box-car, 1: Hanning, 2: Hamming) |
lag |
maximum lag of autocovariance. If |
minmax |
bound for outliers in low side and high side. |
plot |
logical. If |
... |
graphical arguments passed to |
Hanning Window : | = 0.5 |
= 0.25 |
Hamming Window : | = 0.54 |
= 0.23 |
An object of class "spg"
, which is a list with the following components:
period |
periodogram(or raw spectrum). |
smoothed.period |
smoothed log-periodogram. Smoothed periodogram is given if there is a negative value in the smoothed periodogram. |
log.scale |
if |
tsname |
the name of the univariate time series |
Kitagawa, G. (2020) Introduction to Time Series Modeling with Applications in R. Chapman & Hall/CRC.
## BLSALLFOOD data data(BLSALLFOOD) period(BLSALLFOOD) ## seismic Data data(MYE1F) # smoothed periodogram period(MYE1F) # periodogram period(MYE1F, window = 0) # raw spectrum period(MYE1F, window = 0, lag = 200) # Hamming window period(MYE1F, window = 2)
## BLSALLFOOD data data(BLSALLFOOD) period(BLSALLFOOD) ## seismic Data data(MYE1F) # smoothed periodogram period(MYE1F) # periodogram period(MYE1F, window = 0) # raw spectrum period(MYE1F, window = 0, lag = 200) # Hamming window period(MYE1F, window = 2)
Trend estimation by particle filter and smoother.
pfilter(y, m = 10000, model = 0, lag = 20, initd = 0, sigma2, tau2, alpha = 0.99, bigtau2 = NULL, init.sigma2 = 1, xrange = NULL, seed = NULL, plot = TRUE, ...)
pfilter(y, m = 10000, model = 0, lag = 20, initd = 0, sigma2, tau2, alpha = 0.99, bigtau2 = NULL, init.sigma2 = 1, xrange = NULL, seed = NULL, plot = TRUE, ...)
y |
univariate time series. |
||||||||||
m |
number of particles. |
||||||||||
model |
model for the system noise.
|
||||||||||
lag |
lag length for fixed-lag smoothing. |
||||||||||
initd |
type of initial state distribution.
|
||||||||||
sigma2 |
observation noise variance |
||||||||||
tau2 |
system noise variance |
||||||||||
alpha |
mixture weight |
||||||||||
bigtau2 |
variance of the second component |
||||||||||
init.sigma2 |
variance for |
||||||||||
xrange |
specify the lower and upper bounds of the distribution's range. |
||||||||||
seed |
arbitrary positive integer to generate a sequence of uniform random numbers. The default seed is based on the current time. |
||||||||||
plot |
logical. If |
||||||||||
... |
graphical arguments passed to the |
This function performs particle filtering and smoothing for the first order trend model;
|
(system model) |
|
(observation model) |
where is a time series,
is the state vector.
The system noise
and the observation noise
are assumed to be white noises
which follow a Gaussian distribution or a Cauchy distribution, and non-Gaussian distribution, respectively.
The algorithm of the particle filter and smoother are presented in Kitagawa (2020). For more details, please refer to Kitagawa (1996) and Doucet et al. (2001).
An object of class "pfilter"
which has a plot
method. This is a
list with the following components:
llkhood |
log-likelihood. |
||||||||||||
smooth.dist |
marginal smoothed distribution of the trend
|
Kitagawa, G. (1996) Monte Carlo filter and smoother for non-Gaussian nonlinear state space models, J. of Comp. and Graph. Statist., 5, 1-25.
Doucet, A., de Freitas, N. and Gordon, N. (2001) Sequential Monte Carlo Methods in Practice, Springer, New York.
Kitagawa, G. (2020) Introduction to Time Series Modeling with Applications in R. Chapman & Hall/CRC.
pfilterNL
performs particle filtering and smoothing for nonlinear
non-Gaussian state-space model.
data(PfilterSample) y <- PfilterSample ## Not run: pfilter(y, m = 100000, model = 0, lag = 20, initd = 0, sigma2 = 1.048, tau2 = 1.4e-2, xrange = c(-4, 4), seed = 2019071117) pfilter(y, m = 100000, model = 1, lag = 20 , initd = 0, sigma2 = 1.045, tau2 = 3.53e-5, xrange = c(-4, 4), seed = 2019071117) pfilter(y, m = 100000, model = 2, lag = 20 , initd = 0, sigma2 = 1.03, tau2 = 0.00013, alpha = 0.991, xrange = c(-4, 4), seed = 2019071117) ## End(Not run)
data(PfilterSample) y <- PfilterSample ## Not run: pfilter(y, m = 100000, model = 0, lag = 20, initd = 0, sigma2 = 1.048, tau2 = 1.4e-2, xrange = c(-4, 4), seed = 2019071117) pfilter(y, m = 100000, model = 1, lag = 20 , initd = 0, sigma2 = 1.045, tau2 = 3.53e-5, xrange = c(-4, 4), seed = 2019071117) pfilter(y, m = 100000, model = 2, lag = 20 , initd = 0, sigma2 = 1.03, tau2 = 0.00013, alpha = 0.991, xrange = c(-4, 4), seed = 2019071117) ## End(Not run)
Trend estimation by particle filter and smoother via nonlinear state-space model.
pfilterNL(y, m = 10000, lag = 20, sigma2, tau2, xrange = NULL, seed = NULL, plot = TRUE, ...)
pfilterNL(y, m = 10000, lag = 20, sigma2, tau2, xrange = NULL, seed = NULL, plot = TRUE, ...)
y |
univariate time series. |
m |
number of particles. |
lag |
lag length for fixed-lag smoothing. |
sigma2 |
observation noise variance. |
tau2 |
system noise variance. |
xrange |
specify the lower and upper bounds of the distribution's range. |
seed |
arbitrary positive integer to generate a sequence of uniform random numbers. The default seed is based on the current time. |
plot |
logical. If |
... |
graphical arguments passed to the |
This function performs particle filtering and smoothing for the following nonlinear state-space model;
|
(system model) |
|
(observation model) |
where is a time series,
is the state vector.
The system noise
and the observation noise
are assumed to be white noises
which follow a Gaussian distribution and
~
.
The algorithm of the particle filtering and smoothing are presented in Kitagawa (2020). For more details, please refer to Kitagawa (1996) and Doucet et al. (2001).
An object of class "pfilter"
which has a plot
method. This is a
list with the following components:
llkhood |
log-likelihood. |
||||||||||||
smooth.dist |
marginal smoothed distribution of the trend
|
Kitagawa, G. (1996) Monte Carlo filter and smoother for non-Gaussian nonlinear state space models, J. of Comp. and Graph. Statist., 5, 1-25.
Doucet, A., de Freitas, N. and Gordon, N. (2001) Sequential Monte Carlo Methods in Practice, Springer, New York.
Kitagawa, G. (2020) Introduction to Time Series Modeling with Applications in R. Chapman & Hall/CRC.
pfilter
performs particle filtering and smoothing for linear
non-Gaussian state-space model.
data(NLmodel) x <- NLmodel[, 2] pfilterNL(x, m = 100000, lag = 20 , sigma2 = 10.0, tau2 = 1.0, xrange = c(-20, 20), seed = 2019071117)
data(NLmodel) x <- NLmodel[, 2] pfilterNL(x, m = 100000, lag = 20 , sigma2 = 10.0, tau2 = 1.0, xrange = c(-20, 20), seed = 2019071117)
An artificially generated sample data with shifting mean value.
data(PfilterSample)
data(PfilterSample)
A time series of 400 observations.
This data generated by the following models;
~ , |
|
= 0, | 1 <= n <= 100 |
= 1, | 101 <= n <= 200 | ||
= -1, | 201 <= n <= 300 | ||
= 0, | 301 <= n <= 400 |
Plot original data and transformed data with minimum AIC.
## S3 method for class 'boxcox' plot(x, rdata = NULL, ...)
## S3 method for class 'boxcox' plot(x, rdata = NULL, ...)
x |
an object of class |
rdata |
original data, if necessary. |
... |
further graphical parameters may also be supplied as arguments. |
Plot original data and fitted trigonometric polynomial returned by
lsqr
.
## S3 method for class 'lsqr' plot(x, rdata = NULL, ...)
## S3 method for class 'lsqr' plot(x, rdata = NULL, ...)
x |
an object of class |
rdata |
original data, if necessary. |
... |
further graphical parameters may also be supplied as arguments. |
Plot the smoothed density function returned by ngsmth
.
## S3 method for class 'ngsmth' plot(x, type = c("trend", "smt"), theta = 0, phi = 15, expand = 1, col = "lightblue", ticktype= "detail", ...)
## S3 method for class 'ngsmth' plot(x, type = c("trend", "smt"), theta = 0, phi = 15, expand = 1, col = "lightblue", ticktype= "detail", ...)
x |
an object of class |
type |
plotted values, either or both of "trend" and "smt". |
theta , phi , expand , col , ticktype
|
graphical parameters in perspective
plot |
... |
further graphical parameters may also be supplied as arguments. |
Plot trend component of fitted polynomial returned by polreg
.
## S3 method for class 'polreg' plot(x, rdata = NULL, ...)
## S3 method for class 'polreg' plot(x, rdata = NULL, ...)
x |
an object of class |
rdata |
original data, if necessary. |
... |
further graphical parameters may also be supplied as arguments. |
Plot trend component, seasonal component, AR component and noise returned by
season
.
## S3 method for class 'season' plot(x, rdata = NULL, ...)
## S3 method for class 'season' plot(x, rdata = NULL, ...)
x |
an object of class |
rdata |
original data, if necessary. |
... |
further graphical parameters may also be supplied as arguments. |
Plot simulated data of Gaussian / non-Gaussian generated by state space model.
## S3 method for class 'simulate' plot(x, use = NULL, ...)
## S3 method for class 'simulate' plot(x, use = NULL, ...)
x |
an object of class |
use |
start and end time |
... |
further graphical parameters may also be supplied as arguments. |
Plot posterior distribution (mean and standard deviations) of the smoother returned by
tsmooth
.
## S3 method for class 'smooth' plot(x, rdata = NULL, ...)
## S3 method for class 'smooth' plot(x, rdata = NULL, ...)
x |
an object of class |
rdata |
original data, if necessary. |
... |
further graphical parameters may also be supplied as arguments. |
Plot smoothed periodogram or logarithm of smoothed periodogram.
## S3 method for class 'spg' plot(x, type = "vl", ...)
## S3 method for class 'spg' plot(x, type = "vl", ...)
x |
|
type |
type of plot. ("l": lines, "vl" : vertical lines) |
... |
further graphical parameters may also be supplied as arguments. |
Plot trend component and residuals returned by trend
.
## S3 method for class 'trend' plot(x, rdata = NULL, ...)
## S3 method for class 'trend' plot(x, rdata = NULL, ...)
x |
an object of class |
rdata |
original data, if necessary. |
... |
further graphical parameters may also be supplied as arguments. |
Plot evolutionary power spectra obtained by time varying AR model returned
by tvspc
.
## S3 method for class 'tvspc' plot(x, tvv = NULL, dx = 2, dy = 0.25, ...)
## S3 method for class 'tvspc' plot(x, tvv = NULL, dx = 2, dy = 0.25, ...)
x |
an object of class |
tvv |
time varying variance as returned by |
dx |
step width for the X axis. |
dy |
step width for the Y axis. |
... |
further graphical parameters may also be supplied as arguments. |
# seismic data data(MYE1F) v <- tvvar(MYE1F, trend.order = 2, tau2.ini = 6.6e-06, delta = 1.0e-06, plot = FALSE ) z <- tvar(v$nordata, trend.order = 2, ar.order = 8, span = 20, outlier = c(630, 1026), tau2.ini = 6.6e-06, delta = 1.0e-06, plot = FALSE) spec <- tvspc(z$arcoef, z$sigma2, span = 20, nf = 400) plot(spec, tvv = v$tvv, dx = 2, dy = 0.10)
# seismic data data(MYE1F) v <- tvvar(MYE1F, trend.order = 2, tau2.ini = 6.6e-06, delta = 1.0e-06, plot = FALSE ) z <- tvar(v$nordata, trend.order = 2, ar.order = 8, span = 20, outlier = c(630, 1026), tau2.ini = 6.6e-06, delta = 1.0e-06, plot = FALSE) spec <- tvspc(z$arcoef, z$sigma2, span = 20, nf = 400) plot(spec, tvv = v$tvv, dx = 2, dy = 0.10)
Estimate the trend using the AIC best polynomial regression model.
polreg(y, order, plot = TRUE, ...)
polreg(y, order, plot = TRUE, ...)
y |
a univariate time series. |
order |
maximum order of polynomial regression. |
plot |
logical. If |
... |
graphical arguments passed to |
An object of class "polreg"
, which is a list with the following
components:
order.maice |
MAICE (minimum AIC estimate) order. |
sigma2 |
residual variance of the model with order |
aic |
AIC of the model with order |
daic |
AIC - minimum AIC. |
coef |
regression coefficients ( |
trend |
trend component. |
Kitagawa, G. (2020) Introduction to Time Series Modeling with Applications in R. Chapman & Hall/CRC.
# The daily maximum temperatures for Tokyo data(Temperature) polreg(Temperature, order = 7) # Wholesale hardware data data(WHARD) y <- log10(WHARD) polreg(y, order = 15)
# The daily maximum temperatures for Tokyo data(Temperature) polreg(Temperature, order = 7) # Wholesale hardware data data(WHARD) y <- log10(WHARD) polreg(y, order = 15)
Number of rainy days in two years (1975-1976) at Tokyo, Japan.
data(Rainfall)
data(Rainfall)
Integer-valued time series of 366 observations.
The data were obtained from Tokyo District Meteorological Observatory. https://www.data.jma.go.jp/obd/stats/etrn/
Seasonal adjustment by state space modeling.
season(y, trend.order = 1, seasonal.order = 1, ar.order = 0, trade = FALSE, period = NULL, tau2.ini = NULL, filter = c(1, length(y)), predict = length(y), arcoef.ini = NULL, log = FALSE, log.base = "e", minmax = c(-1.0e+30, 1.0e+30), plot = TRUE, ...)
season(y, trend.order = 1, seasonal.order = 1, ar.order = 0, trade = FALSE, period = NULL, tau2.ini = NULL, filter = c(1, length(y)), predict = length(y), arcoef.ini = NULL, log = FALSE, log.base = "e", minmax = c(-1.0e+30, 1.0e+30), plot = TRUE, ...)
y |
a univariate time series with or without the tsp attribute. |
||||||||||
trend.order |
trend order (0, 1, 2 or 3). |
||||||||||
seasonal.order |
seasonal order (0, 1 or 2). |
||||||||||
ar.order |
AR order (0, 1, 2, 3, 4 or 5). |
||||||||||
trade |
logical; if |
||||||||||
period |
If the tsp attribute of
|
||||||||||
tau2.ini |
initial estimate of variance of the system noise |
||||||||||
filter |
a numerical vector of the form |
||||||||||
predict |
the end position of prediction ( |
||||||||||
arcoef.ini |
initial estimate of AR coefficients (for |
||||||||||
log |
logical. If |
||||||||||
log.base |
the letter "e" (default) or "10" specifying the base of logarithmic transformation. Valid only if log = TRUE. |
||||||||||
minmax |
lower and upper limits of observations. |
||||||||||
plot |
logical. If |
||||||||||
... |
graphical arguments passed to |
An object of class "season"
, which is a list with the following
components:
tau2 |
variance of the system noise. |
sigma2 |
variance of the observational noise. |
llkhood |
log-likelihood of the model. |
aic |
AIC of the model. |
trend |
trend component (for |
seasonal |
seasonal component (for |
arcoef |
AR coefficients (for |
ar |
AR component (for |
day.effect |
trading day effect (for |
noise |
noise component. |
cov |
covariance matrix of smoother. |
For time series with the tsp attribute, set frequency
to
period
.
Kitagawa, G. (2020) Introduction to Time Series Modeling with Applications in R. Chapman & Hall/CRC.
# BLSALLFOOD data data(BLSALLFOOD) season(BLSALLFOOD, trend.order = 2, seasonal.order = 1, ar.order = 2) season(BLSALLFOOD, trend.order = 2, seasonal.order = 1, ar.order = 2, filter = c(1, 132)) # Wholesale hardware data data(WHARD) season(WHARD, trend.order = 2, seasonal.order = 1, ar.order = 0, trade = TRUE, log = TRUE) season(WHARD, trend.order = 2, seasonal.order = 1, ar.order = 0, trade = TRUE, filter = c(1, 132), log = TRUE)
# BLSALLFOOD data data(BLSALLFOOD) season(BLSALLFOOD, trend.order = 2, seasonal.order = 1, ar.order = 2) season(BLSALLFOOD, trend.order = 2, seasonal.order = 1, ar.order = 2, filter = c(1, 132)) # Wholesale hardware data data(WHARD) season(WHARD, trend.order = 2, seasonal.order = 1, ar.order = 0, trade = TRUE, log = TRUE) season(WHARD, trend.order = 2, seasonal.order = 1, ar.order = 0, trade = TRUE, filter = c(1, 132), log = TRUE)
Simulate time series by Gaussian State Space Model.
simssm(n = 200, trend = NULL, seasonal.order = 0, seasonal = NULL, arcoef = NULL, ar = NULL, tau1 = NULL, tau2 = NULL, tau3 = NULL, sigma2 = 1.0, seed = NULL, plot = TRUE, ...)
simssm(n = 200, trend = NULL, seasonal.order = 0, seasonal = NULL, arcoef = NULL, ar = NULL, tau1 = NULL, tau2 = NULL, tau3 = NULL, sigma2 = 1.0, seed = NULL, plot = TRUE, ...)
n |
the number of data generated by simulation. |
trend |
initial values of trend component of length |
seasonal.order |
order of seasonal component model (0, 1, 2). |
seasonal |
if |
arcoef |
AR coefficients. |
ar |
initial values of AR component. |
tau1 |
variance of trend component model. |
tau2 |
variance of AR component model. |
tau3 |
variance of seasonal component model. |
sigma2 |
variance of the observation noise. |
seed |
arbitrary positive integer to generate a sequence of uniform random numbers. The default seed is based on the current time. |
plot |
logical. If |
... |
graphical arguments passed to |
An object of class "simulate"
, giving simulated data of Gaussian state
space model.
Kitagawa, G. (2020) Introduction to Time Series Modeling with Applications in R. Chapman & Hall/CRC.
# BLSALLFOOD data data(BLSALLFOOD) m1 <- 2; m2 <- 1; m3 <- 2 z <- season(BLSALLFOOD, trend.order = m1, seasonal.order = m2, ar.order = m3) nl <- length(BLSALLFOOD) trend <- z$trend[m1:1] arcoef <- z$arcoef period <- 12 seasonal <- z$seasonal[(period-1):1] ar <- z$ar[m3:1] tau1 <- z$tau2[1] tau2 <- z$tau2[2] tau3 <- z$tau2[3] simssm(n = nl, trend, seasonal.order = m2, seasonal, arcoef, ar, tau1, tau2, tau3, sigma2 = z$sigma2, seed = 333)
# BLSALLFOOD data data(BLSALLFOOD) m1 <- 2; m2 <- 1; m3 <- 2 z <- season(BLSALLFOOD, trend.order = m1, seasonal.order = m2, ar.order = m3) nl <- length(BLSALLFOOD) trend <- z$trend[m1:1] arcoef <- z$arcoef period <- 12 seasonal <- z$seasonal[(period-1):1] ar <- z$ar[m3:1] tau1 <- z$tau2[1] tau2 <- z$tau2[2] tau3 <- z$tau2[3] simssm(n = nl, trend, seasonal.order = m2, seasonal, arcoef, ar, tau1, tau2, tau3, sigma2 = z$sigma2, seed = 333)
Yearly numbers of sunspots from to 1749 to 1979.
data(Sunspot)
data(Sunspot)
A time series of 231 observations; yearly from 1749 to 1979.
Sunspot is a part of the dataset sunspot.year
from 1700 to 1988.
Value "0" is converted into "0.1" for log transformation.
The daily maximum temperatures in Tokyo (from 1979-01-01 to 1980-04-30).
data(Temperature)
data(Temperature)
A time series of 486 observations.
The data were obtained from Tokyo District Meteorological Observatory. https://www.data.jma.go.jp/obd/stats/etrn/
Estimate the trend by state space model.
trend(y, trend.order = 1, tau2.ini = NULL, delta, plot = TRUE, ...)
trend(y, trend.order = 1, tau2.ini = NULL, delta, plot = TRUE, ...)
y |
a univariate time series. |
trend.order |
trend order. |
tau2.ini |
initial estimate of variance of the system noise |
delta |
search width (for |
plot |
logical. If |
... |
graphical arguments passed to |
The trend model can be represented by a state space model
where ,
and
are matrices with appropriate dimensions.
We assume that
and
are white noises that have
the normal distributions
and
,
respectively.
An object of class "trend"
, which is a list with the following
components:
trend |
trend component. |
residual |
residuals. |
tau2 |
variance of the system noise |
sigma2 |
variance of the observational noise |
llkhood |
log-likelihood of the model. |
aic |
AIC. |
cov |
covariance matrix of smoother. |
Kitagawa, G. (2020) Introduction to Time Series Modeling with Applications in R. Chapman & Hall/CRC.
# The daily maximum temperatures for Tokyo data(Temperature) trend(Temperature, trend.order = 1, tau2.ini = 0.223, delta = 0.001) trend(Temperature, trend.order = 2)
# The daily maximum temperatures for Tokyo data(Temperature) trend(Temperature, trend.order = 1, tau2.ini = 0.223, delta = 0.001) trend(Temperature, trend.order = 2)
Predict and interpolate time series based on state space model by Kalman filter.
tsmooth(y, f, g, h, q, r, x0 = NULL, v0 = NULL, filter.end = NULL, predict.end = NULL, minmax = c(-1.0e+30, 1.0e+30), missed = NULL, np = NULL, plot = TRUE, ...)
tsmooth(y, f, g, h, q, r, x0 = NULL, v0 = NULL, filter.end = NULL, predict.end = NULL, minmax = c(-1.0e+30, 1.0e+30), missed = NULL, np = NULL, plot = TRUE, ...)
y |
a univariate time series |
f |
state transition matrix |
g |
matrix |
h |
matrix |
q |
system noise variance |
r |
observational noise variance |
x0 |
initial state vector |
v0 |
initial state covariance matrix |
filter.end |
end point of filtering. |
predict.end |
end point of prediction. |
minmax |
lower and upper limits of observations. |
missed |
start position of missed intervals. |
np |
number of missed observations. |
plot |
logical. If |
... |
graphical arguments passed to |
The linear Gaussian state space model is
where is a univariate time series,
is an
-dimensional state vector.
,
and
are
,
matrices and a vector of length
, respectively.
is
matrix and
is a scalar.
is system noise and
is observation noise,
where we assume that
,
and
. User should give all the matrices
of a state space model and its parameters. In current version,
,
,
,
,
should be
time invariant.
An object of class "smooth"
, which is a list with the following
components:
mean.smooth |
mean vectors of the smoother. |
cov.smooth |
variance of the smoother. |
esterr |
estimation error. |
llkhood |
log-likelihood. |
aic |
AIC. |
Kitagawa, G. (2020) Introduction to Time Series Modeling with Applications in R. Chapman & Hall/CRC.
Kitagawa, G. and Gersch, W. (1996) Smoothness Priors Analysis of Time Series. Lecture Notes in Statistics, No.116, Springer-Verlag.
## Example of prediction (AR model) data(BLSALLFOOD) BLS120 <- BLSALLFOOD[1:120] z1 <- arfit(BLS120, plot = FALSE) tau2 <- z1$sigma2 # m = maice.order, k=1 m1 <- z1$maice.order arcoef <- z1$arcoef[[m1]] f <- matrix(0.0e0, m1, m1) f[1, ] <- arcoef if (m1 != 1) for (i in 2:m1) f[i, i-1] <- 1 g <- c(1, rep(0.0e0, m1-1)) h <- c(1, rep(0.0e0, m1-1)) q <- tau2[m1+1] r <- 0.0e0 x0 <- rep(0.0e0, m1) v0 <- NULL s1 <- tsmooth(BLS120, f, g, h, q, r, x0, v0, filter.end = 120, predict.end = 156) s1 plot(s1, BLSALLFOOD) ## Example of interpolation of missing values (AR model) z2 <- arfit(BLSALLFOOD, plot = FALSE) tau2 <- z2$sigma2 # m = maice.order, k=1 m2 <- z2$maice.order arcoef <- z2$arcoef[[m2]] f <- matrix(0.0e0, m2, m2) f[1, ] <- arcoef if (m2 != 1) for (i in 2:m2) f[i, i-1] <- 1 g <- c(1, rep(0.0e0, m2-1)) h <- c(1, rep(0.0e0, m2-1)) q <- tau2[m2+1] r <- 0.0e0 x0 <- rep(0.0e0, m2) v0 <- NULL tsmooth(BLSALLFOOD, f, g, h, q, r, x0, v0, missed = c(41, 101), np = c(30, 20))
## Example of prediction (AR model) data(BLSALLFOOD) BLS120 <- BLSALLFOOD[1:120] z1 <- arfit(BLS120, plot = FALSE) tau2 <- z1$sigma2 # m = maice.order, k=1 m1 <- z1$maice.order arcoef <- z1$arcoef[[m1]] f <- matrix(0.0e0, m1, m1) f[1, ] <- arcoef if (m1 != 1) for (i in 2:m1) f[i, i-1] <- 1 g <- c(1, rep(0.0e0, m1-1)) h <- c(1, rep(0.0e0, m1-1)) q <- tau2[m1+1] r <- 0.0e0 x0 <- rep(0.0e0, m1) v0 <- NULL s1 <- tsmooth(BLS120, f, g, h, q, r, x0, v0, filter.end = 120, predict.end = 156) s1 plot(s1, BLSALLFOOD) ## Example of interpolation of missing values (AR model) z2 <- arfit(BLSALLFOOD, plot = FALSE) tau2 <- z2$sigma2 # m = maice.order, k=1 m2 <- z2$maice.order arcoef <- z2$arcoef[[m2]] f <- matrix(0.0e0, m2, m2) f[1, ] <- arcoef if (m2 != 1) for (i in 2:m2) f[i, i-1] <- 1 g <- c(1, rep(0.0e0, m2-1)) h <- c(1, rep(0.0e0, m2-1)) q <- tau2[m2+1] r <- 0.0e0 x0 <- rep(0.0e0, m2) v0 <- NULL tsmooth(BLSALLFOOD, f, g, h, q, r, x0, v0, missed = c(41, 101), np = c(30, 20))
Estimate time varying coefficients AR model.
tvar(y, trend.order = 2, ar.order = 2, span, outlier = NULL, tau2.ini = NULL, delta, plot = TRUE)
tvar(y, trend.order = 2, ar.order = 2, span, outlier = NULL, tau2.ini = NULL, delta, plot = TRUE)
y |
a univariate time series. |
trend.order |
trend order (1 or 2). |
ar.order |
AR order. |
span |
local stationary span. |
outlier |
positions of outliers. |
tau2.ini |
initial estimate of variance of the system noise |
delta |
search width. |
plot |
logical. If |
The time-varying coefficients AR model is given by
where is
-lag AR coefficient at time
and
is a zero mean white noise.
The time-varying spectrum can be plotted using AR coefficient arcoef
and variance of the observational noise sigma2
by tvspc
.
arcoef |
time varying AR coefficients. |
sigma2 |
variance of the observational noise |
tau2 |
variance of the system noise |
llkhood |
log-likelihood of the model. |
aic |
AIC. |
parcor |
PARCOR. |
Kitagawa, G. (2020) Introduction to Time Series Modeling with Applications in R. Chapman & Hall/CRC.
Kitagawa, G. and Gersch, W. (1996) Smoothness Priors Analysis of Time Series. Lecture Notes in Statistics, No.116, Springer-Verlag.
Kitagawa, G. and Gersch, W. (1985) A smoothness priors time varying AR coefficient modeling of nonstationary time series. IEEE trans. on Automatic Control, AC-30, 48-56.
# seismic data data(MYE1F) z <- tvar(MYE1F, trend.order = 2, ar.order = 8, span = 20, outlier = c(630, 1026), tau2.ini = 6.6e-06, delta = 1.0e-06) z spec <- tvspc(z$arcoef, z$sigma2) plot(spec)
# seismic data data(MYE1F) z <- tvar(MYE1F, trend.order = 2, ar.order = 8, span = 20, outlier = c(630, 1026), tau2.ini = 6.6e-06, delta = 1.0e-06) z spec <- tvspc(z$arcoef, z$sigma2) plot(spec)
Estimate evolutionary power spectra by time varying AR model.
tvspc(arcoef, sigma2, var = NULL, span = 20, nf = 200)
tvspc(arcoef, sigma2, var = NULL, span = 20, nf = 200)
arcoef |
time varying AR coefficients. |
sigma2 |
variance of the observational noise. |
var |
time varying variance. |
span |
local stationary span. |
nf |
number of frequencies in evaluating power spectrum. |
return an object of class "tvspc"
giving power spectra, which has a
plot
method (plot.tvspc
).
Kitagawa, G. (2020) Introduction to Time Series Modeling with Applications in R. Chapman & Hall/CRC.
Kitagawa, G. and Gersch, W. (1996) Smoothness Priors Analysis of Time Series. Lecture Notes in Statistics, No.116, Springer-Verlag.
Kitagawa, G. and Gersch, W. (1985) A smoothness priors time varying AR coefficient modeling of nonstationary time series. IEEE trans. on Automatic Control, AC-30, 48-56.
# seismic data data(MYE1F) z <- tvar(MYE1F, trend.order = 2, ar.order = 8, span = 20, outlier = c(630, 1026), tau2.ini = 6.6e-06, delta = 1.0e-06) spec <- tvspc(z$arcoef, z$sigma2) plot(spec)
# seismic data data(MYE1F) z <- tvar(MYE1F, trend.order = 2, ar.order = 8, span = 20, outlier = c(630, 1026), tau2.ini = 6.6e-06, delta = 1.0e-06) spec <- tvspc(z$arcoef, z$sigma2) plot(spec)
Estimate time-varying variance.
tvvar(y, trend.order, tau2.ini = NULL, delta, plot = TRUE, ...)
tvvar(y, trend.order, tau2.ini = NULL, delta, plot = TRUE, ...)
y |
a univariate time series. |
trend.order |
trend order. |
tau2.ini |
initial estimate of variance of the system noise |
delta |
search width. |
plot |
logical. If |
... |
graphical arguments passed to the |
Assuming that , we define a transformed time series
by
where is a Gaussian white noise with mean
and variance
.
is distributed as a
distribution with
degrees of freedom, so the probability density function of
is given by
By further transformation
the probability density function of is given by
Therefore, the transformed time series is given by
where is a double exponential distribution with probability density
function
In the space state model
by identifying trend components of , the log variance of original
time series
is obtained.
An object of class "tvvar"
which has a plot
method. This is a
list with the following components:
tvv |
time varying variance. |
nordata |
normalized data. |
sm |
transformed data. |
trend |
trend. |
noise |
residuals. |
tau2 |
variance of the system noise. |
sigma2 |
variance of the observational noise. |
llkhood |
log-likelihood of the model. |
aic |
AIC. |
tsname |
the name of the univariate time series |
Kitagawa, G. (2020) Introduction to Time Series Modeling with Applications in R. Chapman & Hall/CRC.
Kitagawa, G. and Gersch, W. (1996) Smoothness Priors Analysis of Time Series. Lecture Notes in Statistics, No.116, Springer-Verlag.
Kitagawa, G. and Gersch, W. (1985) A smoothness priors time varying AR coefficient modeling of nonstationary time series. IEEE trans. on Automatic Control, AC-30, 48-56.
# seismic data data(MYE1F) tvvar(MYE1F, trend.order = 2, tau2.ini = 6.6e-06, delta = 1.0e-06)
# seismic data data(MYE1F) tvvar(MYE1F, trend.order = 2, tau2.ini = 6.6e-06, delta = 1.0e-06)
Compute autocovariance and autocorrelation function of the univariate time series.
unicor(y, lag = NULL, minmax = c(-1.0e+30, 1.0e+30), plot = TRUE, ...)
unicor(y, lag = NULL, minmax = c(-1.0e+30, 1.0e+30), plot = TRUE, ...)
y |
a univariate time series. |
lag |
maximum lag. Default is |
minmax |
thresholds for outliers in low side and high side. |
plot |
logical. If |
... |
graphical arguments passed to the |
An object of class "unicor"
which has a plot
method. This is a
list with the following components:
acov |
autocovariances. |
acor |
autocorrelations. |
acov.err |
error bound for autocovariances. |
acor.err |
error bound for autocorrelations. |
mean |
mean of |
tsname |
the name of the univariate time series |
Kitagawa, G. (2020) Introduction to Time Series Modeling with Applications in R. Chapman & Hall/CRC.
# Yaw rate, rolling, pitching and rudder angle of a ship data(HAKUSAN) Yawrate <- HAKUSAN[, 1] unicor(Yawrate, lag = 50) # seismic data data(MYE1F) unicor(MYE1F, lag = 50)
# Yaw rate, rolling, pitching and rudder angle of a ship data(HAKUSAN) Yawrate <- HAKUSAN[, 1] unicor(Yawrate, lag = 50) # seismic data data(MYE1F) unicor(MYE1F, lag = 50)
The monthly record of wholesale hardware data. (January 1967 - November 1979)
data(WHARD)
data(WHARD)
A time series of 155 observations.
The data were obtained from the United States Bureau of Labor Statistics (BLS).