In the first part I created the data for testing the Astronomical/Astrological Hypotheses. In this part, I started by fitting a simple linear regression model.

mod.lm = lm(div.a.ts~vulcan.ts)
summary(mod.lm)
Call:
lm(formula = div.a.ts ~ vulcan.ts)

Residuals:
    Min      1Q  Median      3Q     Max 
-159.30  -53.88  -10.37   53.48  194.05 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 621.23955   20.24209  30.690   <2e-16 ***
vulcan.ts     0.06049    0.05440   1.112    0.273    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 77.91 on 38 degrees of freedom
Multiple R-squared:  0.03151,	Adjusted R-squared:  0.006026 
F-statistic: 1.236 on 1 and 38 DF,  p-value: 0.2731

Clearly, there is little reason to think that the jump points are related to the adjusted divorce counts given the simple regression model. However, since the data is time series data, there is a possibility that the errors are autocorrelated.

The Box-Jenkins Approach

Following the Box-Jenkins approach to fitting time series, I decided to start by looking at the residuals from the adjusted divorce count model as a stationary time series. I generated two plots; a plot of the correlations between present and lagged variables using the R function acf() and a plot of the partial correlations between present and lagged variables using the R function pacf(). For the autocorrelation and the partial autocorrelation plots the blue lines represent the 95% level.

acf(mod.lm$residuals, 
main="autocorrelation of the residuals fitting divorce counts on Ceres, Juno, Vesta, and Pallas breaks")

pacf(mod.lm$residuals, 
main="partial autocorrelation of the residuals
fitting divorce counts on Ceres, Juno, Vesta, and Pallas breaks
")

Here there are the plot: First the auto-correlations and next the partial auto-correlations.

Following the Box-Jenkins approach, the correlation plot will indicate if the time series is a moving average. If there are spikes in the correlation plot, the spikes indicate the orders of the moving average. An autoregressive time series will have an exponentially decaying correlation plot.

The partial correlation plot indicates the order of the process if the process is autoregressive. The order of the auto regressive process is given by the value of the lag just before the partial correlation goes to essentially zero.

A process with both moving average and auto regressive terms is hard to identify using correlation and partial correlation plots, but the process is indicated by decay that is not exponential, but starts after a few lags (see the Wikipedia page Box-Jenkins).

Looking at the above plots, the autocorrelation plot shows negative spikes at lag 1 and lag 13 and is large in absolute value at lags 3, 12, and 14, but otherwise is consistent with random noise. The partial correlation plot has negative spikes at lags 1, 2, and 13 and is large in absolute value at lag 15, indicating that the process may not be stationary. Neither plot decays.

The Periodogram and Normalized Cumulative Periodogram

Using periodograms is another way you can evaluate autocorrelation. A periodogram can find periodicities in a time series. Below are plots of the periodogram and the normalized cumulative periodogram created using spec.pgram() and cpgram() in R. In the periodogram, there are spikes at about the frequencies .35 and .40 that are quite large. The frequencies correspond to lags of 2.9 and 2.5, somewhat close to the lags of 2 and 3 months found in the correlation plot. The normalized cumulative periodogram shows a departure from random noise at a frequency of around .29, or a lag of around 3.4. The blue lines in the normalized cumulative periodogram plot are the 95% confidence limits for random noise.

Plot the periodogram and cumulative periodogram of the residuals from the model.

spec.pgram(mod.lm$residuals, 
main="Residuals from Simple Regression
Raw Periodogram")

cpgram(mod.lm$residuals, 
main="Residuals from Simple Regression
Cumulative Periodogram")

First the periodogram and next the cumulative periodogram.

Note also that at the frequency of twenty, corresponding to a lag of two, the periodogram is also high.

Using arima() in R to Compare Models

I used the function arima() in R to fit some autoregressive moving average models with orders up to thirteen, both with and without the astronomical regressor. The best model, using the Aktaike information criteria (see the Wikipedia page Aktaike information criterion), was $$Y{_t}=\beta_{0} + \beta_{1}X_{t} + \theta{_1}\epsilon_{t-1} + \theta{_3} \epsilon_{t-1} + \theta_{13} \epsilon_{t-13} + \epsilon_{t}$$ where Y_t is the number of divorces in month t, X_t is the number of days since the last shift in the mean of the celestial longitudes of Ceres, Pallas, Vesta, and Juno over the shortest arc, and ε_t is the difference between the expected number of divorces and the actual number of divorces for month t.

The estimated coefficients, with the respective standard errors, are given in the table below.

The final model was fit using arima(). Below, two models are fit, the final one (with the regressor) and one without the regressor.

In arima(), in the argument fixed, an NA indicates the parameter is fit and a 0 indicates no parameter is fit. In fixed there is a value for every value indicated by the arguments order and xreg, plus the intercept. So in the first model below, there are 13 moving average terms, three of which are fit, and there is an intercept and one regressor term, for 15 values in fixed.

The model with the regressor.

mod.arima.r.1 = arima(div.a.ts, 
                 order=c(0,0,13), 
                 fixed=c(NA,0,NA,rep(0,9),
                 NA,NA,NA), 
                 xreg=vulcan.ts) 
mod.arima.r.1
Call:
arima(x = div.a.ts, order = c(0, 0, 13), xreg = vulcan.ts, fixed = c(NA, 0, 
    NA, rep(0, 9), NA, NA, NA))

Coefficients:
          ma1  ma2     ma3  ma4  ma5  ma6  ma7  ma8  ma9  ma10  ma11  ma12     ma13  intercept
      -0.6822    0  0.3855    0    0    0    0    0    0     0     0     0  -0.5975   623.7144
s.e.   0.2744    0  0.1799    0    0    0    0    0    0     0     0     0   0.2270    10.4553
      vulcan.ts
         0.0541
s.e.     0.0309

sigma^2 estimated as 2618:  log likelihood = -220.54,  aic = 453.07

The model without the regressor

mod.arima.r.2 = arima(div.a.ts, 
                 order=c(0,0,13), 
                 fixed=c(NA,0,NA,rep(0,9),
                 NA,NA))
mod.arima.r.2
Call:
arima(x = div.a.ts, order = c(0, 0, 13), fixed = c(NA, 0, NA, rep(0, 9), NA, 
    NA))

Coefficients:
          ma1  ma2     ma3  ma4  ma5  ma6  ma7  ma8  ma9  ma10  ma11  ma12     ma13  intercept
      -0.6044    0  0.4800    0    0    0    0    0    0     0     0     0  -0.6596   640.3849
s.e.   0.2025    0  0.1775    0    0    0    0    0    0     0     0     0   0.2266     4.7133

sigma^2 estimated as 2655:  log likelihood = -221.92,  aic = 453.85

The coefficient for the slope of the regression is not at least 1.96 times the size of the standard errors of the coefficient, however, the Aktaike information criteria is the smallest for the models at which I looked. The model without the coefficient is a close contender, with the Aktaike information for the two models being 453.07 and 453.85.

The residuals from the two models are plotted together using ts.plot().

ts.plot(cbind(mod.arima.r.1$residuals, 
     mod.arima.r.2$residuals), col=4:3,
     xlab="Time", 
     ylab="Residuals",
     main=
     "Residuals for Two Moving Average Models")
legend("topright", 
     c("With Regression Coefficient", 
     "Without Regression Coefficient"), 
     fill=4:3, cex=.7)

Below is a plot of the residuals for the two competing models. To me the residuals look more even for the model given above.

The estimates indicate that the intercept is about 624 divorces; that the number of divorces increases by about 5.4 for every one hundred day increase in the time after the midpoint over the shortest arc of the longitudes of Juno, Ceres, Vesta, and Pallas shifts by 90 or 180 degrees; that errors tend to reverse in direction over one lag and thirteen lags, but tend to be in the same direction over three lags.

Some Concluding Diagnostic Plots

Below are some diagnostic plots; the plot of the residuals versus the fitted values, the Q-Q plot of the residuals versus normal variates, and a time series plot of the actual values and the fitted values.

First the graphic screen is set up to plot three plots in a row. Next the residuals are plotted against the fitted values. Next the residuals are plotted against standard normal percentiles. Last the data is plotted with the fitted values using ts.plot().

par(mfrow=c(1,3))
plot(as.numeric(div.a.ts-
    mod.arima.r.1$residuals), 
  as.numeric(mod.arima.r.1$residuals),
  main="Plot of Residuals versus Fitted Values for the Adjusted Divorce Data", 
  xlab="Fitted Value", 
  ylab="Residual", 
  font.main=1)
qqnorm(mod.arima.r.1$residuals, 
  main="QQ Normal Plot for Residuals from the Divorce Model", 
  font.main=1)
ts.plot(div.a.ts,div.a.ts-
    mod.arima.r.1$residuals,
  ylab="Counts", 
  xlab="Year",
  main="Adjusted Divorce Counts and Fitted Counts", 
  col=4:3)
legend("bottomleft", 
  c("Adjusted Divorce Counts","Fitted Adjusted Divorce Counts"), 
  col=4:3, lwd=2, cex=.7)
cor.test(div.a.ts-mod.arima.r.1$residuals, 
  mod.arima.r.1$residuals)

Here it is the plot:

The residual versus fitted value plot looks like there is some correlation between the residuals and the fitted values. Testing the correlation gives a p-value of 0.0112, so the model is not ideal. While the model does not pick up most of the extreme data points, the residual and QQ normal plots do not show any great deviations from a model of normal errors with constant variance. The fitted values were found by subtracting the residuals from the actual data points.

The correlation between the fitted values and the residuals is found and tested for equality with zero using cor.test()

cor.test(div.a.ts-mod.arima.r.1$residuals, 
         mod.arima.r.1$residuals)
Pearson's product-moment correlation

data:  div.a.ts - mod.arima.r.1$residuals and mod.arima.r.1$residuals
t = 2.6651, df = 38, p-value = 0.01123
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.09736844 0.63041834
sample estimates:
      cor 
0.3968411 

The Function to Find the Best Model

Here is the function I created to compare models. The function is set up to only do moving averages. It takes a while to run with ‘or’ equal to 13 for moving averages since 2 to the power of 13, which equals 8,192, arima()’s are run, since each element of fixed can take on either ‘0’ or ‘NA’. If both autoregressive and moving average models are run, then for ‘or’ equal to 13, the function would have to fit 4 to the power of 13 models to get all combinations or 67,108,864 models, which is beyond the capacity of R, at least on my machine.

In the function, ‘ts’ is the time series to be modeled, ‘phi’ is the vector (matrix) of regressor(s), and ‘or’ is the highest moving average order to be included. To run just autoregressive models, in the calls to arima(), move ‘or’ from the third element of ‘order’ to the first element of ‘order’. In order to run both autoregressive and moving average models, remove the pound signs and add pound signs below where the pound signs were.

In the function, first ‘phi’ is made into a matrix, if a vector. Next the matrices ‘fx.mat1’ and ‘fx.mat2’ are created with the vectors to be assigned to ‘fixed’ in the rows. ‘fx.mat1’ contains regressors and ‘fx.mat2’ does not. Next ‘res1’ and ‘res2’ are created to hold the values of the Aktaike criteria for each model; ‘res1’ for models with a regressor and ‘res2’ for models without a regressor. Next the arima models are fit in a for() loop. Then the numeric order of the elements in ‘res1’ and ‘res2’ are found using order() and only the ten models with the smallest Aktaike criteria are kept, for each of ‘res1’ and ‘res2’. The variables ‘st1’ and ‘st2’ contain the indices of the models with the smallest Aktaike criteria. Last, the values for the ten smallest Aktaike criteria and the values for ‘fixed’ for the ten smallest Aktaike criteria are printed, for models both with and without the regressor.

arima.comp=function(ts=div.a.ts, 
phi=vulcan.ts, or=3) {
  
  phi <- as.matrix(phi)
  if (nrow(phi)==1) phi <- t(phi)
  fx.lst <- list(0:1)
  #for (i in 2:(2*or)) fx.lst <- 
  #c(fx.lst,list(0:1))
  for (i in 2:or) fx.lst <- 
    c(fx.lst,list(0:1))
  
  fx.mat <- 
    as.matrix(expand.grid(fx.lst))
  #fx.phi1 <- matrix(NA, 2^(2*or), 
  #  1+ncol(phi))
  fx.phi1 <- matrix(NA, 2^or, 
    1+ncol(phi))
  
  fx.mat[fx.mat==1] <- NA
  fx.mat1 <- 
    cbind(fx.mat,fx.phi1)
  fx.mat2 <- 
    cbind(fx.mat,fx.phi1[,1])
  
  #res1 <- numeric(2^(2*or))
  #res2 <- numeric(2^(2*or))
  
  res1 <- numeric(2^or)
  res2 <- numeric(2^or)
  
  #for (i in 1:2^(2*or)) {
  for (i in 1:2^or) {
    #res1[i] <- arima(ts, order=c(or,0,or), 
    #  xreg=phi, fixed=fx.mat1[i,])$aic
    #res2[i] <- arima(ts, order=c(or,0,or), 
    #  fixed=fx.mat2[i,])$aic
    
    res1[i] <- arima(ts, order=c(0,0,or), 
      xreg=phi, fixed=fx.mat1[i,])$aic
    res2[i] <- arima(ts, order=c(0,0,or), 
      fixed=fx.mat2[i,])$aic
  }
  
  st1 <- order(res1)
  st1 <- st1[1:min(length(st1),10)]
  print(res1[st1])
  print(fx.mat1[st1,])
  
  st2 <- order(res2)
  st2 <- st2[1:min(length(st2),10)]
  print(res2[st2])
  print(fx.mat2[st2,])
}

In the next part, we use the covariance structure of the model selected in this part to estimate standard errors for differing linear combinations of the monthly counts. The results are interpreted from an astrological point of view.

1. Regression model with auto correlated errors – Part 1, the data
2. Linear Regression from Scratch in R
3. R for Publication by Page Piccinini: Lesson 5 – Analysis of Variance (ANOVA)
4. R for Publication by Page Piccinini: Lesson 4 – Multiple Regression
5. R for Publication by Page Piccinini: Lesson 3 – Logistic Regression