Time Series Analysis with State Space Model

Description

R functions for statistical analysis, modeling and simulation of time series with state space model.

Details

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.

References

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).

BLSALLFOOD Data

Description

The monthly time series of the number of workers engaged in food industries in the United States (January 1967 - December 1979).

Usage

data(BLSALLFOOD)

Format

A time series of 156 observations.

Source

The data were obtained from the United States Bureau of Labor Statistics (BLS).

Ship's Navigation Data

Description

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.

Usage

data(HAKUSAN)

Format

A data frame with 1000 observations on the following 4 variables.

Source

The data were offered by Prof. K. Ohtsu of Tokyo University of Marine Science and Technology.

Examples

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)

Haibara Data

Description

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.

Usage

data(Haibara)

Format

A data frame with 400 observations on the following 2 variables.

Source

The data were offered by Dr. M. Takahashi and Dr. N. Matsumoto of National Institute of Advanced Industrial Science and Technology.

Examples

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)

Seismic Data

Description

The time series of East-West components of seismic waves, recorded every 0.02 seconds.

Usage

data(MYE1F)

Format

A time series of 2600 observations.

Source

Takanami, T. (1991), "ISM data 43-3-01: Seismograms of foreshocks of 1982 Urakawa-Oki earthquake", Ann. Inst. Statist. Math., 43, 605.

The Nonlinear State-Space Model Data

Description

The series generated by the nonlinear state-space model.

Usage

data(NLmodel)

Format

A matrix with 100 rows and 2 columns.

Details

The system model x_n and the observation model y_n are generated by following state-space model:

x_n = \frac{1}{2} x_{n-1} + \frac{25x_{n-1}}{x_{n-1}^2+1} + 8cos(1.2n) + v_n

y_n = \frac{x_n^2}{10} + w_n,

where v_n ~ N(0, 1), w_n ~ N(0, 10), v_0 ~ N(0, 5).

Nikkei225

Description

A daily closing values of the Japanese stock price index, Nikkei225, quoted from January 4, 1988, to December 30, 1993.

Usage

data(Nikkei225)

Format

A time series of 1480 observations.

Source

https://indexes.nikkei.co.jp/nkave/archives/data

Sample Data for Particle Filter and Smoother

Description

An artificially generated sample data with shifting mean value.

Usage

data(PfilterSample)

Format

A time series of 400 observations.

Details

This data generated by the following models;

Rainfall Data

Description

Number of rainy days in two years (1975-1976) at Tokyo, Japan.

Usage

data(Rainfall)

Format

Integer-valued time series of 366 observations.

Source

The data were obtained from Tokyo District Meteorological Observatory. https://www.data.jma.go.jp/obd/stats/etrn/

Sunspot Number Data

Description

Yearly numbers of sunspots from to 1749 to 1979.

Usage

data(Sunspot)

Format

A time series of 231 observations; yearly from 1749 to 1979.

Details

Sunspot is a part of the dataset sunspot.year from 1700 to 1988. Value "0" is converted into "0.1" for log transformation.

Temperatures Data

Description

The daily maximum temperatures in Tokyo (from 1979-01-01 to 1980-04-30).

Usage

data(Temperature)

Format

A time series of 486 observations.

Source

The data were obtained from Tokyo District Meteorological Observatory. https://www.data.jma.go.jp/obd/stats/etrn/

Wholesale Hardware Data

Description

The monthly record of wholesale hardware data. (January 1967 - November 1979)

Usage

data(WHARD)

Format

A time series of 155 observations.

Source

The data were obtained from the United States Bureau of Labor Statistics (BLS).

Univariate AR Model Fitting

Description

Fit a univariate AR model by the Yule-Walker method, the least squares (Householder) method or the PARCOR method.

Usage

arfit(y, lag = NULL, method = 1, plot = TRUE, ...)

Arguments

Value

An object of class "arfit" which has a plot method. This is a list with the following components:

References

Kitagawa, G. (2020) Introduction to Time Series Modeling with Applications in R. Chapman & Hall/CRC.

Examples

# Sunspot number data
data(Sunspot)
arfit(log10(Sunspot), lag = 20, method = 1)

# BLSALLFOOD data
data(BLSALLFOOD)
arfit(BLSALLFOOD)

Calculate Characteristics of Scalar ARMA Model

Description

Calculate impulse response function, autocovariance function, autocorrelation function and characteristic roots of given scalar ARMA model.

Usage

armachar(arcoef = NULL, macoef = NULL, v, lag = 50, nf = 200, plot = TRUE, ...)

Arguments

Details

The ARMA model is given by

y_t - a_1y_{t-1} - \dots - a_py_{t-p} = u_t - b_1u_{t-1} - \dots - b_qu_{t-q},

where p is AR order, q is MA order and u_t is a zero mean white noise.

Characteristic roots of AR / MA operator is a list with the following components:

• re: real part R

• im: imaginary part I

• amp: \sqrt{R^2+I^2}

• atan: \arctan(I/R)

• degree

Value

An object of class "arma" which has a plot method. This is a list with components:

References

Kitagawa, G. (2020) Introduction to Time Series Modeling with Applications in R. Chapman & Hall/CRC.

Examples

# 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)

Scalar ARMA Model Fitting

Description

Fit a scalar ARMA model by maximum likelihood method.

Usage

armafit(y, ar.order, ar = NULL, ma.order, ma = NULL)

Arguments

Value

References

Kitagawa, G. (2020) Introduction to Time Series Modeling with Applications in R. Chapman & Hall/CRC.

Examples

# 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)

Scalar ARMA Model Fitting

Description

Estimate all ARMA models within the user-specified maximum order by maximum likelihood method.

Usage

armafit2(y, ar.order, ma.order)

Arguments

Value

References

Kitagawa, G. (2020) Introduction to Time Series Modeling with Applications in R. Chapman & Hall/CRC.

Examples

# Sunspot number data
data(Sunspot)
y <- log10(Sunspot)
armafit2(y, ar.order = 5, ma.order = 5)

Box-Cox Transformation

Description

Compute Box-Cox transformation and find an optimal lambda with minimum AIC.

Usage

boxcox(y, plot = TRUE, ...)

Arguments

Value

An object of class "boxcox", which is a list with the following components:

References

Kitagawa, G. (2020) Introduction to Time Series Modeling with Applications in R. Chapman & Hall/CRC.

Examples

# Sunspot number data
data(Sunspot)
boxcox(Sunspot)

# Wholesale hardware data
data(WHARD)
boxcox(WHARD)

Cross-Covariance and Cross-Correlation

Description

Compute cross-covariance and cross-correlation functions of the multivariate time series.

Usage

crscor(y, lag = NULL, outmin = NULL, outmax = NULL, plot = TRUE, ...)

Arguments

Value

An object of class "crscor" which has a plot method. This is a list with the following components:

References

Kitagawa, G. (2020) Introduction to Time Series Modeling with Applications in R. Chapman & Hall/CRC.

Examples

# 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 via FFT

Description

Compute a periodogram of the univariate time series via FFT.

Usage

fftper(y, window = 1, plot = TRUE, ...)

Arguments

Details

Value

An object of class "spg", which is a list with the following components:

Note

We assume that the length N of the input time series y is a power of 2. If N is not a power of 2, calculate using the FFT by appending 0's behind the data y.

References

Kitagawa, G. (2020) Introduction to Time Series Modeling with Applications in R. Chapman & Hall/CRC.

Examples

# Yaw rate, rolling, pitching and rudder angle of a ship
data(HAKUSAN)
YawRate <- HAKUSAN[, 1]
fftper(YawRate, window = 0)

Kullback-Leibler Information

Description

Compute Kullback-Leibler information.

Usage

klinfo(distg = 1, paramg = c(0, 1), distf = 1, paramf, xmax = 10)

Arguments

Value

References

Kitagawa, G. (2020) Introduction to Time Series Modeling with Applications in R. Chapman & Hall/CRC.

Examples

# 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)

Decomposition of Time Interval to Stationary Subintervals

Description

Decompose time series to stationary subintervals and estimate local spectrum.

Usage

lsar(y, max.arorder = 20, ns0, plot = TRUE, ...)

Arguments

Value

An object of class "lsar" which has a plot method. This is a list with the following components:

References

Kitagawa, G. (2020) Introduction to Time Series Modeling with Applications in R. Chapman & Hall/CRC.

Examples

# seismic data
data(MYE1F)
lsar(MYE1F, max.arorder = 10, ns0 = 100)

Estimation of the Change Point

Description

Precisely estimate a change point of subinterval for locally stationary AR model.

Usage

lsar.chgpt(y, max.arorder = 20, subinterval, candidate, plot = TRUE, ...)

Arguments

Value

An object of class "chgpt" which has a plot method. This is a list with the following components:

References

Kitagawa, G. (2020) Introduction to Time Series Modeling with Applications in R. Chapman & Hall/CRC.

Examples

# 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))

The Least Squares Method via Householder Transformation

Description

Compute regression coefficients of the model with minimum AIC by the least squares method via Householder transformation.

Usage

lsqr(y, lag = NULL, period = 365, plot = TRUE, ...)

Arguments

Value

An object of class "lsqr", which is a list with the following components:

References

Kitagawa, G. (2020) Introduction to Time Series Modeling with Applications in R. Chapman & Hall/CRC.

Examples

# The daily maximum temperatures in Tokyo
data(Temperature)
lsqr(Temperature, lag = 10)

Yule-Walker Method of Fitting Multivariate AR Model

Description

Fit a multivariate AR model by the Yule-Walker method.

Usage

marfit(y, lag = NULL)

Arguments

Value

An object of class "maryule", which is a list with the following components:

References

Kitagawa, G. (2020) Introduction to Time Series Modeling with Applications in R. Chapman & Hall/CRC.

Examples

# 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)

Least Squares Method for Multivariate AR Model

Description

Fit a multivariate AR model by least squares method.

Usage

marlsq(y, lag = NULL)

Arguments

Value

An object of class "marlsq", which is a list with the following components:

References

Kitagawa, G. (2020) Introduction to Time Series Modeling with Applications in R. Chapman & Hall/CRC.

Examples

# 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)

Cross Spectra and Power Contribution

Description

Compute cross spectra, coherency and power contribution.

Usage

marspc(arcoef, v, plot = TRUE, ...)

Arguments

Value

An object of class "marspc" which has a plot method. This is a list with the following components:

References

Kitagawa, G. (2020) Introduction to Time Series Modeling with Applications in R. Chapman & Hall/CRC.

Examples

# 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)

Simulation by Non-Gaussian State Space Model

Description

Simulation by non-Gaussian state space model.

Usage

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, ...)

Arguments

Value

An object of class "simulate", giving simulated data of non-Gaussian state space model.

References

Kitagawa, G. (2020) Introduction to Time Series Modeling with Applications in R. Chapman & Hall/CRC.

Examples

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))

Non-Gaussian Smoothing

Description

Trend estimation by non-Gaussian smoothing.

Usage

ngsmth(y, noisev = 2, tau2, bv = 1.0, noisew = 1, sigma2, bw = 1.0, 
       initd = 1, k = 200, plot = TRUE, ...)

Arguments

Details

Consider a one-dimensional state space model

x_n = x_{n-1} + v_n,

y_n = x_n + w_n,

where the observation noise w_n is assumed to be Gaussian distributed and the system noise v_n is assumed to be distributed as the Pearson system

q(v_n) = c/(\tau^2 + v_n^2)^b

with \frac{1}{2} < b < \infty and c = \tau^{2b-1}\Gamma(b) \bigm/ \Gamma(\frac{1}{2})\Gamma(b-\frac{1}{2}).

This broad family of distributions includes the Cauchy distribution (b = 1) and t-distribution (b = (k+1)/2).

Value

An object of class "ngsmth", which is a list with the following components:

References

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.

Examples

## 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)

Probability Density Function

Description

Evaluate probability density function for normal distribution, Cauchy distribution, Pearson distribution, exponential distribution, Chi-square distributions, double exponential distribution and uniform distribution.

Usage

pdfunc(model = "norm", mean = 0, sigma2 = 1, mu = 0, tau2 = 1, shape,
       lambda = 1, side = 1, df, xmin = 0, xmax = 1, plot = TRUE, ...)

Arguments

Value

An object of class "pdfunc" which has a plot method. This is a list with the following components:

References

Kitagawa, G. (2020) Introduction to Time Series Modeling with Applications in R. Chapman & Hall/CRC.

Examples

# 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

Description

Compute a periodogram of the univariate time series.

Usage

period(y, window = 1, lag = NULL, minmax = c(-1.0e+30, 1.0e+30),
       plot = TRUE, ...)

Arguments

Details

Value

An object of class "spg", which is a list with the following components:

References

Kitagawa, G. (2020) Introduction to Time Series Modeling with Applications in R. Chapman & Hall/CRC.

Examples

## 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)

Particle Filtering and Smoothing

Description

Trend estimation by particle filter and smoother.

Usage

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, ...)

Arguments

Details

This function performs particle filtering and smoothing for the first order trend model;

where y_n is a time series, x_n is the state vector. The system noise v_n and the observation noise w_n 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).

Value

An object of class "pfilter" which has a plot method. This is a list with the following components:

References

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.

See Also

pfilterNL performs particle filtering and smoothing for nonlinear non-Gaussian state-space model.

Examples

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)

Particle Filtering and Smoothing for Nonlinear State-Space Model

Description

Trend estimation by particle filter and smoother via nonlinear state-space model.

Usage

pfilterNL(y, m = 10000, lag = 20, sigma2, tau2, xrange = NULL, seed = NULL,
          plot = TRUE, ...)

Arguments

Details

This function performs particle filtering and smoothing for the following nonlinear state-space model;

where y_n is a time series, x_n is the state vector. The system noise v_n and the observation noise w_n are assumed to be white noises which follow a Gaussian distribution and v_0 ~ N(0, 5).

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).

Value

An object of class "pfilter" which has a plot method. This is a list with the following components:

References

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.

See Also

pfilter performs particle filtering and smoothing for linear non-Gaussian state-space model.

Examples

data(NLmodel)
x <- NLmodel[, 2]
pfilterNL(x, m = 100000, lag = 20 , sigma2 = 10.0, tau2 = 1.0,
          xrange = c(-20, 20), seed = 2019071117)

Plot Box-Cox Transformed Data

Description

Plot original data and transformed data with minimum AIC.

Usage

## S3 method for class 'boxcox'
plot(x, rdata = NULL, ...)

Arguments

Plot Fitted Trigonometric Polynomial

Description

Plot original data and fitted trigonometric polynomial returned by lsqr.

Usage

## S3 method for class 'lsqr'
plot(x, rdata = NULL, ...)

Arguments

Plot Smoothed Density Function

Description

Plot the smoothed density function returned by ngsmth.

Usage

## S3 method for class 'ngsmth'
plot(x, type = c("trend", "smt"), theta = 0, phi = 15,
          expand = 1, col = "lightblue", ticktype= "detail", ...)

Arguments

Plot Fitted Polynomial Trend

Description

Plot trend component of fitted polynomial returned by polreg.

Usage

## S3 method for class 'polreg'
plot(x, rdata = NULL, ...)

Arguments

Plot Trend, Seasonal and AR Components

Description

Plot trend component, seasonal component, AR component and noise returned by season.

Usage

## S3 method for class 'season'
plot(x, rdata = NULL, ...)

Arguments

Plot Simulated Data Generated by State Space Model

Description

Plot simulated data of Gaussian / non-Gaussian generated by state space model.

Usage

## S3 method for class 'simulate'
plot(x, use = NULL, ...)

Arguments

Plot Posterior Distribution of Smoother

Description

Plot posterior distribution (mean and standard deviations) of the smoother returned by tsmooth.

Usage

## S3 method for class 'smooth'
plot(x, rdata = NULL, ...)

Arguments

Plot Smoothed Periodogram

Description

Plot smoothed periodogram or logarithm of smoothed periodogram.

Usage

## S3 method for class 'spg'
plot(x, type = "vl", ...)

Arguments

Plot Trend and Residuals

Description

Plot trend component and residuals returned by trend.

Usage

## S3 method for class 'trend'
plot(x, rdata = NULL, ...)

Arguments

Plot Evolutionary Power Spectra Obtained by Time Varying AR Model

Description

Plot evolutionary power spectra obtained by time varying AR model returned by tvspc.

Usage

## S3 method for class 'tvspc'
plot(x, tvv = NULL, dx = 2, dy = 0.25, ...)

Arguments

Examples

# 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)

Polynomial Regression Model

Description

Estimate the trend using the AIC best polynomial regression model.

Usage

polreg(y, order, plot = TRUE, ...)

Arguments

Value

An object of class "polreg", which is a list with the following components:

References

Kitagawa, G. (2020) Introduction to Time Series Modeling with Applications in R. Chapman & Hall/CRC.

Examples

# The daily maximum temperatures for Tokyo
data(Temperature)
polreg(Temperature, order = 7)

# Wholesale hardware data
data(WHARD)
y <- log10(WHARD)
polreg(y, order = 15)

Seasonal Adjustment

Description

Seasonal adjustment by state space modeling.

Usage

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, ...)

Arguments

Value

An object of class "season", which is a list with the following components:

Note

For time series with the tsp attribute, set frequency to period.

References

Kitagawa, G. (2020) Introduction to Time Series Modeling with Applications in R. Chapman & Hall/CRC.

Examples

# 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)

Simulation by Gaussian State Space Model

Description

Simulate time series by Gaussian State Space Model.

Usage

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, ...)

Arguments

Value

An object of class "simulate", giving simulated data of Gaussian state space model.

References

Kitagawa, G. (2020) Introduction to Time Series Modeling with Applications in R. Chapman & Hall/CRC.

Examples

# 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)

Trend Estimation

Description

Estimate the trend by state space model.

Usage

trend(y, trend.order = 1, tau2.ini = NULL, delta, plot = TRUE, ...)

Arguments

Details

The trend model can be represented by a state space model

x_n = Fx_{n-1} + Gv_n,

y_n = Hx_n + w_n,

where F, G and H are matrices with appropriate dimensions. We assume that v_n and w_n are white noises that have the normal distributions N(0,\tau^2) and N(0, \sigma^2), respectively.

Value

An object of class "trend", which is a list with the following components:

References

Kitagawa, G. (2020) Introduction to Time Series Modeling with Applications in R. Chapman & Hall/CRC.

Examples

# 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)

Prediction and Interpolation of Time Series

Description

Predict and interpolate time series based on state space model by Kalman filter.

Usage

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, ...)

Arguments

Details

The linear Gaussian state space model is

x_n = F_n x_{n-1} + G_n v_n,

y_n = H_n x_n + w_n,

where y_n is a univariate time series, x_n is an m-dimensional state vector.

F_n, G_n and H_n are m \times m, m \times k matrices and a vector of length m , respectively. Q_n is k \times k matrix and R_n is a scalar. v_n is system noise and w_n is observation noise, where we assume that E(v_n, w_n) = 0, v_n \sim N(0, Q_n) and w_n \sim N(0, R_n). User should give all the matrices of a state space model and its parameters. In current version, F_n, G_n, H_n, Q_n, R_n should be time invariant.

Value

An object of class "smooth", which is a list with the following components:

References

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.

Examples

## 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))

Time Varying Coefficients AR Model

Description

Estimate time varying coefficients AR model.

Usage

tvar(y, trend.order = 2, ar.order = 2, span, outlier = NULL, tau2.ini = NULL,
     delta, plot = TRUE)

Arguments

Details

The time-varying coefficients AR model is given by

y_t = a_{1,t}y_{t-1} + \ldots + a_{p,t}y_{t-p} + u_t

where a_{i,t} is i-lag AR coefficient at time t and u_t 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.

Value

References

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.

See Also

tvspc, plot.tvspc

Examples

# 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)

Evolutionary Power Spectra by Time Varying AR Model

Description

Estimate evolutionary power spectra by time varying AR model.

Usage

tvspc(arcoef, sigma2, var = NULL, span = 20, nf = 200)

Arguments

Value

return an object of class "tvspc" giving power spectra, which has a plot method (plot.tvspc).

References

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.

Examples

# 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)

Time Varying Variance

Description

Estimate time-varying variance.

Usage

tvvar(y, trend.order, tau2.ini = NULL, delta, plot = TRUE, ...)

Arguments

Details

Assuming that \sigma_{2m-1}^2 = \sigma_{2m}^2, we define a transformed time series s_1,\dots,s_{N/2} by

s_m = y_{2m-1}^2 + y_{2m}^2,

where y_n is a Gaussian white noise with mean 0 and variance \sigma_n^2. s_m is distributed as a \chi^2 distribution with 2 degrees of freedom, so the probability density function of s_m is given by

f(s) = \frac{1}{2\sigma^2} e^{-s/2\sigma^2}.

By further transformation

z_m = \log \left( \frac{s_m}{2} \right),

the probability density function of z_m is given by

g(z) = \frac{1}{\sigma^2} \exp{ \left\{ z-\frac{e^z}{\sigma^2} \right\} } = \exp{ \left\{ (z-\log\sigma^2) - e^{(z-\log\sigma^2)} \right\} }.

Therefore, the transformed time series is given by

z_m = \log \sigma^2 + w_m,

where w_m is a double exponential distribution with probability density function

h(w) = \exp{\{w-e^w\}}.

In the space state model

z_m = t_m + w_m

by identifying trend components of z_m, the log variance of original time series y_n is obtained.

Value

An object of class "tvvar" which has a plot method. This is a list with the following components:

References

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.

Examples

# seismic data
data(MYE1F)
tvvar(MYE1F, trend.order = 2, tau2.ini = 6.6e-06, delta = 1.0e-06)

Autocovariance and Autocorrelation

Description

Compute autocovariance and autocorrelation function of the univariate time series.

Usage

unicor(y, lag = NULL, minmax = c(-1.0e+30, 1.0e+30), plot = TRUE, ...)

Arguments

Value

An object of class "unicor" which has a plot method. This is a list with the following components:

References

Kitagawa, G. (2020) Introduction to Time Series Modeling with Applications in R. Chapman & Hall/CRC.

Examples

# 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)