In jlivsey/sigex: Ecce Signum: Multivariate Signal Extraction

knitr::opts_knit$set(root.dir = 'C:/Users/neide/Documents/GitHub/sigex')

Non-Defense Capitalization

This illustration examines monthly Shipments and New Orders from the Manufacturing, Shipments, and Inventories survey. (This is seasonally adjusted monthly data covering the sample period January 1992 to April 2020, downloaded July 21, 2020 (4:45 PM), U.S. Census Bureau, obtained from https://www.census.gov/mtis/index.html by selecting Non-Defense Capital Goods, and either Value of Shipments or New Orders.) Call this ndc for short. The data for New Orders is not available at January 1992, since this series starts at February 1992; so this value is entered as an NA. This is an example of ragged edge data.

Loading Packages

This code loads the data into a variable named ndc.

## wipe
rm(list=ls())
library(devtools)
library(Rcpp)

load_all(".")
root.dir <- getwd()

Enter Metadata

• We perform a bivariate analysis by fitting a VAR model.
• First enter metadata, including starting date and frequency.

start.date <- c(1992,1)
end.date <- c(2020,5)
period <- 12

• Next enter names of the series.
• The last argument of sigex.load generates a time series plot.

dataALL.ts <- sigex.load(ndc,start.date,period,
                         c("Shipments","NewOrders"),TRUE)

• On the basis of this plot, we decide (below) to utilize no transformation for the data.

Select Spans and Transforms

• We have the choice of either log or no log.
• The aggregate option will sum across the indicated series, in case we want to analyze an aggregate. Here we set it to FALSE.
• We can also select a subcomponent of series with subseries.
• The range argument of sigex.prep will select a subset of dates for the time series to be analyzed.

transform <- "none"
aggregate <- FALSE
subseries <- c(1,2)
begin.date <- start(dataALL.ts)
end.date <- end(dataALL.ts)
range <- NULL
data.ts <- sigex.prep(dataALL.ts,transform,aggregate,subseries,range,TRUE)

Model Declaration

• Now we want to define our VAR model.
• We begin by defining dimension $N$ and sample size $T$.

N <- dim(data.ts)[2]
T <- dim(data.ts)[1]

Preliminary Construction

• To define the VAR model, we begin with fitting via Yule-Walker, using base R.
• To do this, we are careful to omit the first observation (due to the missing value), and difference the time series to remove trend non-stationarity.

ar.fit <- ar.yw(diff(ts(ndc[2:T,])))
p.order <- ar.fit$order

• The parameter estimates will be used as initial values for our MLE routine.
• But we must first arrange the VAR parameter estimates in the right order, which requires some permuting.
• Also the innovation covariance matrix is broken into its Cholesky decomposition using getGCD. This function outputs a list object, including the lower triangular entry of $L$ and the log of the diagonal entries of $D$.

par.yw <- aperm(ar.fit$ar,c(2,3,1))
covmat.yw <- getGCD(ar.fit$var.pred,2)

• Next, the array of VAR parameters are mapped via var.par2pre to a pre-parameter vector var.out of length 28.
• Then psi.init is assembled from the Yule-Walker fit, including the sample means, which are incorporated as the last two pre-parameters.

var.out <- var2.par2pre(par.yw)
psi.init <- as.vector(c(covmat.yw[[1]][2,1],log(covmat.yw[[2]]),
                        var.out,colMeans(diff(ts(ndc[2:T,])))))

Model Construction

• Now we define the model.
• The list object mdl stores all our specifications about the model, but none of the parameter estimates.

mdl <- NULL
mdl <- sigex.add(mdl,seq(1,N),"varma",c(p.order,0),NULL,"process",c(1,-1))

• The call to sigex.add indicates there is one main component.
• The second argument indicates the rank configuration for the process; setting vrank equal to seq(1,N) ensures a full rank covariance matrix.
• The third argument gives the model type as varma;
  the order $p=7$, $q=0$ of the VARMA is given in the fourth argument.
• The fifth argument provides parameter bounds, and does not apply to this type of model, and hence is set to NULL.
• The sixth argument associates a label.
• The last argument is for the scalar differencing operator $\delta^{(1)} (z)$, expressed as a vector of coefficients.

mdl <- sigex.meaninit(mdl,data.ts,0)

• There are no regressors specified, but we still need to call sigex.meaninit so that the time trend polynomial is set up correctly.
• The last argument gives the order of this time trend polynomial. The value $0$ indicates there is a constant mean effect assumed.

Model Estimation

• Now we will fit the model, and check the fit.
• First we need to initialize. The variable psi.mle correspond to the pre-parameters, and get mapped to par.mle, the parameters, via sigex.psi2par.

constraint <- NULL
psi.mle <- psi.init
par.mle <- sigex.psi2par(psi.mle,mdl,data.ts)

MLE Fitting

• Maximum Likelihood Estimation is done using BFGS. We use the divergence, which is $-2$ times the log Gaussian likelihood, with constants removed.
• With the debug option set to TRUE, the values of the divergence are printed to the console.

fit.mle <- sigex.mlefit(data.ts,par.mle,constraint,mdl,"bfgs",debug=FALSE)
psi.mle <- sigex.eta2psi(fit.mle[[1]]$par,constraint)
hess <- fit.mle[[1]]$hessian
par.mle <- fit.mle[[2]]

• This yields a divergence of r sigex.lik(psi.mle,mdl,data.ts).
• The MLE pre-parameter is given next.

print(round(psi.mle,digits=4))

Residual Analysis

• We obtain, format, and store the residuals.
• Then we compute the autocorrelations, and plot.

resid.mle <- sigex.resid(psi.mle,mdl,data.ts)[[1]]
resid.mle <- sigex.load(t(resid.mle),start(data.ts),frequency(data.ts),
                        colnames(data.ts),TRUE)
resid.acf <- acf(resid.mle,lag.max=4*period,plot=TRUE)$acf

• We can also examine the condition numbers, to see if there is rank reduction possible.
• Typically there won't be any rank reduction with a single component model.

log(sigex.conditions(data.ts,psi.mle,mdl))

Model Checking

• We can examine the Portmanteau statistic for residual serial correlation. We use $48$ lags.

sigex.portmanteau(resid.mle,4*period,length(psi.mle))

• The output indicates there is little residual autocorrelation, which is confirmed by inspection of the sample autocorrelation.
• We can also inspect normality via the Shapiro-Wilks normality test.

sigex.gausscheck(resid.mle)

• There is some non-normality indicated. However, we are not principally interested in the marginal distributions of our time series model.
• We can also examine the Hessian matrix, and compute the t statistics for pre-parameters.

print(eigen(hess)$values)

• If the Hessian is not positive definite (this can happen if nonlinear optimization terminates wrongfully at a saddlepoint), then the standard error is set to zero.
• In this case there is no problem.

tstats <- sigex.tstats(mdl,psi.mle,hess,constraint)
print(tstats)

• These results could be used to insert zeroes in the pre-parameters, although they need not correspond to zeroes in the parameters.

Bundle

• Having completed our model analysis, we store our results in a single list object.

analysis.mle <- sigex.bundle(data.ts,transform,mdl,psi.mle)

Casting

• Now the modeling and fitting is completed, we can do forecasting and midcasting.
• First we reload from the bundle.
• If we had other fitted models in another bundle, we could load that instead.

data.ts <- analysis.mle[[1]]
mdl <- analysis.mle[[3]]
psi <- analysis.mle[[4]]
param <- sigex.psi2par(psi,mdl,data.ts)

• Next, we generate $50$ aftcasts and forecasts.
• Midcast for the missing value at the first time point is automatically generated.
• Uncertainty is also computed.

window <- 50
data.casts <- sigex.midcast(psi,mdl,data.ts,window)
extract.casts <- sigex.castextract(data.ts,data.casts,mdl,window,param)

• The output of sigex.midcast is a list with two items: the matrix of casts, and a 4-array of casting error covariances.
• Note that these casts are for the series with all fixed effects removed.
• This output can be formatted for use via sigex.castextract, which inserts the casts in the correct places in the time series and adds an extended trend regression effect (it is extended forward and backward by the castspan value).
• The output of sigex.castextract can be utilized by the graphics routine.

castcol <- "black"
fade <- 60
dataPad.ts <- rbind(matrix(NA,nrow=window,ncol=N),data.ts,matrix(NA,nrow=window,ncol=N))
par(mfrow=c(2,1))
for(i in 1:N)
{
  plot(ts(dataPad.ts[,i],start=start.date,frequency=period),
       xlab="Year",ylab="",lwd=1,col=1,
       ylim=c(min(extract.casts[[3]][,i]),max(extract.casts[[2]][,i])))
  sigex.graph(extract.casts,NULL,start.date,period,i,0,castcol,fade)
}
dev.off()

• This code sets up grey shading for cast uncertainty.
• For each component time series, the raw data is plotted after first padding with NA before and after.

jlivsey/sigex documentation built on March 20, 2024, 3:17 a.m.