Economic impact of presidents and war

Spencer Graves and Jouni Helske

2016-05-01

Abstract

This vignette illustrates the use of the KFAS package for R to evaluate the claim that World War II (WW II), not Franklin Roosevelt (FDR), ended the Great Depression. This is the dominant view of the evolution of the US economy between 1933 and 1945. It is, however, superficially contradicted by our naive review of a plot of average annual income in the US (real US Gross Domestic Product (GDP) per capita, 1790-2014, from MeasuringWorth) in two ways: First, other wars have not been accompanied by such obvious growth spurts. Second, FDR’s economic performance even without World War II seems unmatched by any other U.S. president. This document describes alternative attempts to quantify the economic impact of wars separate from that of the administration (“president”). This includes a detailed discussion and comparison of alternative models of economic growth and how KFAS was used to fit them. We have three conclusions: The models we fit failed to find a signficant macroeconomic impact of war but did find that the performance of the U.S. economy under the Hoover and FDR administrations was dramatically different from that for any other presidency. Second, these conclusions must be considered tentative, because none of these models considered other variables such as inflation and unemployment. Third, more generally, KFAS provides a powerful set of tools for modeling extensions to the present analysis, e.g., mixing data from multiple sources such as joint models of GDP per capita with unemployment and mixing annual and quarterly data.

Introduction

The Legacy of the New Deal and the Great Depression still carry profound implications for the current economic policies of the world’s most economically advanced countries. Horowitz (2011) recently claimed that Hoover had more aggressive deficit spending and public works programs than Franklin Roosevelt, and far from contributing to the recovery, these programs actually made the economy worse. He concluded that “The persistence of the Hoover myth continues to justify the counterproductive policies of the Obama administration and thereby prevents markets from generating the economic recovery of which they are fully capable.”

McChesney (2014) said, “In earlier eras it had been assumed that there was an economic ‘guns and butter’ tradeoff, and that military spending had to occur at the expense of other sectors of the economy. However, one of the lessons of the economic expansion in Nazi Germany, followed by the experience of the United States itself in arming for the Second World War, was that big increases in military spending could act as huge stimulants to the economy.” Pollin and Garrett-Peltier (2009) report four alternatives that would produce more jobs than military spending: Tax cuts for household consumption, clean energy, health care and education would produce 27, 47, 69 and 151 percent more jobs, respectively, than military spending.

Beyond this, the work of Milton Friedman and others suggests that modeling like this should also consider inflation and unemployment.

Whatever the specific cause, the growth of the U.S. economy under FDR before the war was not matched by that of any other president (not counting the negative record of Herbert Hoover). This document describes alternative attempts to quantify the economic impact of wars separate from that of the administration (“president”). This includes a detailed discussion and comparison of alternative models of economic growth and how KFAS was used to fit them.

The best of the models we fit include a “Presidency” effect. Models with a “war” effect failed to fit the data any better than without any such effect. However, these conclusions must be considered tentative, pending further research using other data and possibly other ways of modeling the macroeconomic impact of war.

Initial Plots

We start by plotting the MeasuringWorth estimates of average annual income in the US (adjusted for inflation) from 1790 to 2014:

library(Ecdat)
## Loading required package: Ecfun
## 
## Attaching package: 'Ecdat'
## The following object is masked from 'package:datasets':
## 
##     Orange
plot(realGDPperCapita/1000~Year, USGDPpresidents)
real US GDP per capita (thousands of 2009 $)

real US GDP per capita (thousands of 2009 $)

This plot could be changed in many ways to make it more informative. First, this data.frame covers many years for which we have no GDP data. Let’s start by checking missing data:

which(GDPna <- is.na(USGDPpresidents$realGDPperCapita))
##  [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23
## [24] 24 25 26 27 28 29 30 31 32 33

The only NAs are the first 33 observations. Those contain data on other variables for years prior to 1790. We won’t use any of that data for this analysis. Let’s discard it:

GDP <- USGDPpresidents[!GDPna,]

Let’s try that plot again, this time on a log scale with improved labeling:

plot(realGDPperCapita/1000~Year, GDP, log='y', type='l',
     ylab='thousands of 2009 dollars', las=1, 
     main='Real US GDP per capita')
real US GDP per capita (thousands of 2009 $)

real US GDP per capita (thousands of 2009 $)

This plot suggests relatively steady growth from 1790 to 2014 with the exception of a dramatic rise roughly two thirds of the way through the series preceded and followed by declines.

The president is identified in an ordered factor, “executive”, giving the names of the chief executive (president since 1790). We need to drop the unused levels from the colonial period:

GDP$executive <- ordered(GDP$executive)
levels(GDP$executive)
##  [1] "Washington"     "JAdams"         "Jefferson"      "Madison"       
##  [5] "Monroe"         "JQAdams"        "Jackson"        "VanBuren"      
##  [9] "Tyler"          "Polk"           "Taylor-Filmore" "Pierce"        
## [13] "Buchanan"       "Lincoln"        "AJohnson"       "Grant"         
## [17] "Hayes"          "Arthur"         "Cleveland"      "BHarrison"     
## [21] "Cleveland2"     "McKinley"       "TRoosevelt"     "Taft"          
## [25] "Wilson"         "Harding"        "Coolidge"       "Hoover"        
## [29] "FRoosevelt"     "Truman"         "Eisenhower"     "Kennedy"       
## [33] "LJohnson"       "Nixon"          "Ford"           "Carter"        
## [37] "Reagan"         "GHWBush"        "Clinton"        "GWBush"        
## [41] "Obama"

This data.frame includes three variables for war. The first is an ordered factor. We need to drop the unused levels of “war”, as we just did for “executive”:

GDP$war <- ordered(GDP$war)
levels(GDP$war)
##  [1] ""                            "Northwest Indian War"       
##  [3] "War of 1812"                 "Mexican-American"           
##  [5] "U.S. Civil War"              "Spanish-American-Philippine"
##  [7] "World War I"                 "World War II"               
##  [9] "Korean"                      "Vietnam"

The next war variable is the number of “battleDeaths” for each war apportioned by the number of days in each year between the official start and end of each war. However, it hardly makes sense to use this directly as a measure of the impact of the war on the economy, as the nation grew from 3.9 million souls in 1790 to 319 million in 2014. Instead, we use battle deaths per million population:

plot(battleDeathsPMP~Year, GDP, type='l',
     ylab='battle deaths per million pop', las=1, 
     main='Battle deaths during war')
Battle deaths per million population

Battle deaths per million population

Wars involving fewer than 10 battle deaths per million US population per year are not included. No war since Vietnam has reached that threshold. This explains why the 2001-2014 War in Afghanistan and the 2003-2011 Iraq war do not appear here: As noted in the help page for USGDPpresidents, the two together averaged fewer than 3 battle deaths per year per million population.

This document first adds annotations to the above plot of realGDPperCapita then considers a series of models that estimate a base growth rate separate from other potential explanatory variables such as the administration (the president) and wars. Readers interested in the results are encouraged to focus on the plots and tables, skipping the math and software.

Graphical analysis

We want a simple computation that will identify the exceptional period in the above plot. First, let’s compute the range:

head(GDP, 1)
##    Year  CPI GDPdeflator population.K realGDPperCapita  executive
## 34 1790 8.86        4.34         3929          1107.32 Washington
##                     war battleDeaths battleDeathsPMP Keynes
## 34 Northwest Indian War       115.85           29.49      0
tail(GDP, 1)
##     Year    CPI GDPdeflator population.K realGDPperCapita executive war
## 258 2014 236.74      108.32       319173            50397     Obama    
##     battleDeaths battleDeathsPMP Keynes
## 258            0               0      0
range(GDP$realGDPperCapita)
## [1]  1107.32 50397.00

It’s common to work with log(dollars) rather than dollars directly. This is obvious to many who work with finance but foreign to many others. For the latter group, it’s worth noting that a $500 drop in annual income would have been catastrophic in 1790, representing almost half the average income at the time. However a similar $500 drop around 2014 would be hardly noticeable, being just under 1 percent of the 2014 average income. Moreover, a small change in logarithms can be interpreted as percentage, since log(1+x) is approximately x when x is small (using natural logarithms).

We start looking for this exceptional period by computing the differences in log(real GDP per capita):

difGDP <- diff(log(GDP$realGDPperCapita))

A run up or down ends where difGDP changes sign from one year to the next:

chgSgn <- c(TRUE, (head(difGDP, -1)*tail(difGDP, -1))<0)

We convert this logical vector to an index numbering the runs using cumsum:

GDPrun <- cumsum(chgSgn)

We use tapply to identify the beginning and end of each run and to compute run length and total and average growth:

yr0 <- head(GDP$Year, -1)
runStart <- tapply(yr0, GDPrun, min)
yr1 <- tail(GDP$Year, -1)
runEnd <- tapply(yr1, GDPrun, max)
N <- nrow(GDP)
runLength <- tapply(rep(1, N-1), GDPrun, sum)
totGrowth <- tapply(difGDP, GDPrun, sum)
avgGrowth <- tapply(difGDP, GDPrun, mean)

We organize these numbers into a data.frame for convenience:

str(GDPruns <- data.frame(runStart=runStart, 
      runEnd=runEnd, runLength=runLength, 
      totalGrowth=totGrowth, meanGrowth=avgGrowth))
## 'data.frame':    81 obs. of  5 variables:
##  $ runStart   : int [1:81(1d)] 1790 1796 1797 1802 1803 1806 1808 1814 1817 1818 ...
##   ..- attr(*, "dimnames")=List of 1
##   .. ..$ : chr  "1" "2" "3" "4" ...
##  $ runEnd     : int [1:81(1d)] 1796 1797 1802 1803 1806 1808 1814 1817 1818 1819 ...
##   ..- attr(*, "dimnames")=List of 1
##   .. ..$ : chr  "1" "2" "3" "4" ...
##  $ runLength  : num [1:81(1d)] 6 1 5 1 3 2 6 3 1 1 ...
##   ..- attr(*, "dimnames")=List of 1
##   .. ..$ : chr  "1" "2" "3" "4" ...
##  $ totalGrowth: num [1:81(1d)] 0.2434 -0.0101 0.0934 -0.0141 0.0423 ...
##   ..- attr(*, "dimnames")=List of 1
##   .. ..$ : chr  "1" "2" "3" "4" ...
##  $ meanGrowth : num [1:81(1d)] 0.0406 -0.0101 0.0187 -0.0141 0.0141 ...
##   ..- attr(*, "dimnames")=List of 1
##   .. ..$ : chr  "1" "2" "3" "4" ...

To focus on the potentially most interesting periods, we select the runs with average growth exceeding 0.08 (positive or negative):

GDPruns[abs(GDPruns$meanGrowth)>0.08,]
##    runStart runEnd runLength totalGrowth  meanGrowth
## 39     1894   1895         1  0.09141688  0.09141688
## 44     1907   1908         1 -0.13379206 -0.13379206
## 48     1913   1914         1 -0.09886640 -0.09886640
## 56     1929   1933         4 -0.33649531 -0.08412383
## 57     1933   1937         4  0.33349161  0.08337290
## 59     1938   1944         6  0.66854652  0.11142442

This gives us six runs, three of which occur in the period that jumps off the above plot at us: 1929-1933, 1933-1937, and 1938-1944. The first period is the presidency of Herbert Hoover, and the last two encompass the FDR years of the Great Depression and World War II. (The break in the middle is the “Recession of 1937-38”, which convinced FDR to restore the high levels of government spending that seemed to make the difference between the disaster of 1929-33 and the growth of 1934-36. If we lower the threshold to 0.06, we get 10 runs, the first of which is 1861-1863 at the beginning of the American Civil War. The last is the post WW II recession. The others do not appear to be related to presidents or wars.)

We label these points on the above plot:

plot(realGDPperCapita/1000~Year, GDP, type='l', log='y', 
     main=paste0('US average annual income\n', 
                 '(GDP per capita)'), 
     ylab='thousands of 2009 dollars', las=1)     
abline(v=c(1929, 1933, 1945), lty='dashed', col='grey')
text(1930, 2.5, "Hoover", srt=90, cex=0.9)
text(1939.5, 30, 'FDR', srt=90, cex=1.1, col='blue')
Presidents & exceptional change in GDP

Presidents & exceptional change in GDP

We also want to label the wars on the plot. Since we’ll want war labels on other plots, and since we used multiple lines of code to do it, we’ll create a function to support this:

plotPresWars <- function(x, ..., 
     yHoover=0.8, yFDR=2, yWar=seq(1.7, by=-.1, len=9), 
     cex.war=0.67){
  plot(x, ...)
  abline(v=c(1929, 1933, 1945), lty='dashed')
  usr <- par('usr')
  Hy <- (usr[3]+yHoover*diff(usr[3:4]))
  if(par('ylog'))Hy <- exp(Hy)
  text(1930, Hy, 'Hoover', srt=90, cex=.9, xpd=NA)
  Ry <- (usr[3]+yFDR*diff(usr[3:4]))
  if(par('ylog'))Ry <- exp(Ry)
  text(1939, Ry, 'FDR', srt=90, cex=1.1, col='blue', xpd=NA)
  yrMid <- mean(usr[1:2])
  dyr <- diff(usr[1:2])
  y. <- seq(usr[3], usr[4], len=5)
  dy <- diff(y.[c(2,4)])
  y2.5 <- mean(y.[c(2,3)])
  wars <- levels(GDP$war)
  nWars <- length(wars)
  for(i in 2:nWars){
    w <- wars[i]
    sel <- (GDP$war==w)
    yrs <- range(GDP$Year[sel])
    abline(v=yrs, lty='dotted', col='grey')
    yr. <- mean(yrs)
    wy <- (y2.5 + yWar[i-1]*dy)
    if(par('ylog'))wy <- exp(wy)
    text(yr., wy, w, srt=90, col='red', cex=cex.war, 
         xpd=NA)
  }
  invisible("done")
}

plotPresWars(GDP$realGDPperCapita/1000, log='y', 
             main=paste0('US average annual income\n', 
                         '(GDP per capita)'), 
             ylab='thousands of 2009 dollars', las=1)
Function to plot presidents and wars

Function to plot presidents and wars

This plot has problems. Most importantly, the horizontal axis is labeled “Index” not “Year”. We can fix that by adding “Year” as the first argument to the call to plot. A less obvious ways is to convert “realGDPperCapita” to an object of class ts:

GDP$realGDPperCapita <- ts(GDP$realGDPperCapita, 
                             GDP$Year[1])
plotPresWars(GDP$realGDPperCapita/1000, log='y', 
             main=paste0('US average annual income\n', 
                         '(GDP per capita)'), 
             ylab='thousands of 2009 dollars', las=1)
Presidents, wars & change in GDP

Presidents, wars & change in GDP

The rest of this vignette describes various models of economic growth that include especially the impact of changes of administration and wars. First we outline our basic approach to modeling these kinds of data.

Bayesian sequential updating: overview

The models used here are all state space models, which we discuss in general here, giving specific examples later. State space modeling is also known as Kalman (Kalman filtering, Kalman smoothing), and dynamic (linear) modeling. For present purposes, we use the KFAS package, though many other R packages provide some support for this class of models.

The theory behind these models can be described in terms of a 3-step Bayesian updating cycle to add new observations to a probability distribution summarizing the available knowledge about the state of the system under study (Graves, 2007):

  1. Receive and store new data.

  2. Create a prior distribution for the new data by modifying the posterior computed after the last observation to allow for possible changes in state in the time between the last and current observations.

  3. Combine the new data with the prior to produce a posterior distribution of the state.

In the present discussion, the system under study is the US economy at various times, which we will summarize in a vector of one or more numbers that represent the state or condition of the economy at a particular point in time. (The term “state” in “state space” refers to this type of representation.) The data will be annual observations on real US GDP per capita (realGDPperCapita).

The first model of this type discussed here considers univariate observations \(y_t\) = log(real GDP per capita) on a hidden univariate state \(\alpha_t\), which represents the portion of the annual number that persists. In other words, the differences between the observations and predictions based on the hidden state are assumed to be uncorrelated across time; any serial correlation in \(y_t\) is assumed to be due to the hidden state, \(\alpha_t\).

The general upward trend in the above plots of real GDP per capita indicate that a model with a univariate state will not fit very well. We use it to fix ideas and to serve as a basis for evaluating the importance of estimating the growth rate at time \(t\) separate from the level – and for evaluating the general utility of this approach to modeling time series.

With univariate observations on a univariate state, the result is almost a traditional exponentially weighted moving average (EWMA), except that the weight on the last observation is adjusted with each observation consistent with the Bayesian sequential updating cycle outlined above. In this case, this weight converges to an asymptote determined by the ratio of the migration and noise variances (Graves, Bisgaard, and Kulahci, 2002). In the more general Kalman filtering terminology, the weight on the last observation (or the weight on the deviation of that observation from prediction) is called the Kalman gain. Bayesian sequential updating produces a varying Kalman gain, though in many applications the Kalman gain converges to a constant. Before fitting any model, we first describe the general KFAS model and show how the Bayesian EWMA described above is a special case. This formulation allows the system to vary over time even to the extent of supporting asynchronous observations on different variables at irregular points in time.

The next generalization beyond an EWMA is double exponential smoothing: It has a two-dimensional state vector with the first dimension being the level and the second being the rate of growth. The plots above suggests that double exponential smoothing might work relatively well for log(realGDPperCapita) except for 1929-1947. To model that period, we consider various ways of modeling presidents and wars.

In the following, we shall consider seven different though related models:

Model state vector components time varying components parameters to estimate
1 EWMA1: Standard Bayesian EWMA 1:level 2:lnQ, lnH
2 EWMA2: double EWMA 2:level, slope 3:lnQ.lvl, lnQ.slope, lnH
3 EWMA2pres1: EWMA2 allowing slope to change only between presidents 2:level, slope Q 3:lnQ.lvl, lnQ.slope, lnH
4 EWMA2pres2: EWMA2 assuming no continuity in slope between presidents 2:level, slope Q, T 3:lnQ.lvl, lnQ.slope, lnH
5 EWMA2pres3: EWMA2 with an adjustment for each president 3:level, slope, Pres Q, T 4:lnQ.lvl, lnQ.slope, lnQ.Pres, lnH
6 EWMA2war0: EWMA2 with fixed effect of war 4:level, slope, war0, war1 T 3:lnQ.lvl, lnQ.slope, lnH
7 EWMA2war1: EWMA2 with random effect of war 4:level, slope, war0, war1 Q, T 6:lnQ.lvl, lnQ.slope, lnQ.war00, lnQ.war11, lnQ.war01, lnH

1. The simplest Kalman Filter: an Exponentially Weighted Moving Average (EWMA)

Consider a univariate noisy measurement, \(y_t\), of a univariate state, \(\alpha_t\):

\[y_t = \alpha_t + \epsilon_t\]

where the state, \(\alpha_t\), follows a random walk:

\[\alpha_{t+1} = \alpha_t + \eta_t.\]

In this model, we assume that observation error, \(\epsilon_t\), and the random migration, \(\eta_t\), are normally distributed with mean 0 and variance, \(H\) and \(Q\), respectively, and the initial state \(\alpha_1\) is normal with mean \(a_1\) and variance \(P_1\). We sometimes write \(a_1\) and \(P_1\) as \(a_{1|0}\) and \(P_{1|0}\) to make it clear that these are the prior mean and variance of \(\alpha_1\) before \(y_1\) is available (using the conditional subscript notation, following, e.g., Harvey (1989)). We may also use \(a_{t|t}\) and \(P_{t|t}\) to denote the posterior mean and variance given \(D_t = \{y_t, y_{t-1}, ...\}\). The prior distribution for \(\alpha_t\) given \(D_{t-1}\) and the posterior for \(\alpha_t\) given \(D_t\) will be computed recursively. (See Durbin and Koopman (2012b) for more detail on this model.)

These equations are called the measurement and state transition equations, respectively.

We will estimate \(H\) and \(Q\) to maximize the likelihood. We will also estimate the series \(\alpha_t\) conditional on these estimates \(H\), \(Q\), \(a_1\) and \(P_1\).

With this notation, we can provide more details in discussing the 3-step Bayesian sequential updating process outlined above:

  1. A new observation, \(y_t\), arrives.

  2. The prior distribution \(f(\alpha_t | D_{t-1})\) is computed from the previous posterior \(f(\alpha_{t-1} | D_{t-1})\), where \(D_{t-1}\) is all the information available prior to the arrival of \(y_t\); for \(t\) = 0, this is already assumed as \(\alpha_1 \tilde\ N(a_{1|0}, P_{1|0})\). Otherwise, we have by recursion that the posterior for \((\alpha_{t-1} | D_{t-1})\) is \(N(a_{t-1|t-1}, P_{t-1|t-1})\). From this using a linear state transition equation, we get that the next prior \((\alpha_t | D_{t-1})\) is \(N(a_{t|t-1}, P_{t|t-1})\). In particular, the state transition equation above gives us \(a_{t|t-1} = a_{t-1|t-1}\) and \(P_{t|t-1} = P_{t-1|t-1} + Q.\) (The double subscripts help clarify the math, which can be confusing otherwise.)

  3. Bayes’ theorem is then used to update this prior to incorporate the information in the new observation. This converts \(f(\alpha_t | D_{t-1})\) into \(f(\alpha_t | D_t)\).

The distribution of the state vector \(\alpha_t\) is governed by certain hyperparameters, suppressed in the notation \(f(\alpha_t | D_t)\). These hyperparameters are estimated by maximum likelihood.

KFAS package for Kalman Filtering and Smoothing

State space modeling with the KFAS package assumes the following generalization of the above measurement and state transition equations:

\[y_t = Z_t \alpha_t + \epsilon_t\]

\[\alpha_{t+1} = T_t \alpha_t + R_t \eta_t.\]

where \(\epsilon_t \tilde\ N(0, H_t)\), \(\eta_t \tilde\ N(0, Q_t)\) and \(\alpha_1 \tilde\ N(a_1, P_1)\), independent of each other. In this, \(y_t\) is a \(p\)-vector (where \(p\) may actually vary with \(t\)), \(\alpha_t\) is \(m\)-dimensional, and \(\eta_t\) has dimension \(k\). (See Durbin and Koopman (2012a) for more detail on this model.)

KFAS modeling begins with a formula that is converted to an “SSModel” object using the SSModel function. This is passed to fitSSM to estimate unknown parameters to maximize logLik.SSModel. Confidence and prediction intervals can be obtained using predict.SSModel. Function KFS will produce filtered and smoothed estimates of the state vector in components “a” and “alphahat”, respectively, of an object of class “KFS”:

formula -> SSModel -> fitSSM -> KFS

Let’s write a suitable formula for our Bayesian EWMA:

library(KFAS)
EWMA1fmla <- log(realGDPperCapita)~
  SSMtrend(1,Q=list(matrix(NA)), a1=log(realGDPperCapita[1]) )

This is converted to an SSModel as follows:

EWMA1mdl <-SSModel(EWMA1fmla, GDP, H=matrix(NA) )

This model has 2 unknowns, the NA’s in EWMA1fmla and SSModel:

The default behavior of fitSSM is to estimate the logarithms of these two variance. fitSSM requires us to provide starting values for log(Q) and log(H), in that order:

EWMA1fit <- fitSSM(EWMA1mdl, inits=c(lnQ=0, lnH=0))

Unfortunately, there are two standard likelihoods commonly used for time series programmed into KFAS: diffuse and marginal. By default, the fitSSM function maximizes the diffuse likelihood, because it’s faster to compute and is optimized with the same parameter values in many cases. However, if \(Z_t\) or \(T_t\) are functions of the unknown parameters, this may not be true; see logLik.SSModel. In such cases, the marginal likelihood should be used.

In the examples considered here, neither \(Z_t\) nor \(T_t\) will vary with unknown parameters. However, we will need the marginal likelihood, because it does not change with two different forumlations of the same model and the diffuse likelihood can change. Thus, with nested models, the marginal likelihood should never be less for the larger model, but the diffuse likelihood might. We will see that later.

For now, let’s check this claim by maximizing the marginal likelihood:

EWMA1fitM <- fitSSM(EWMA1mdl, inits=EWMA1fit$optim.out$par, 
                    marginal=TRUE)
EWMA1fit$optim.out$par
##        lnQ        lnH 
##  -6.167604 -35.447595
EWMA1fitM$optim.out$par
##        lnQ        lnH 
##  -6.167604 -35.447595
logLik(EWMA1fit$model)
## [1] 372.9323
logLik(EWMA1fitM$model, marginal=TRUE)
## [1] 375.6403

As expected.

Helske (2016b) notes that the diffuse log-likelihood can be written as follows:

\[\log(L_d) = -\frac{1}{2}\sum^{n}_{t=1}\sum^{p_t}_{i=1}w_{i,t} ,\]

where

and \(F_{\infty,i,t}>0\), \(v^2_{i,t}\), and \(F_{i,t}\) are components of the list output by fitSSM. Thus, we can manually compute the diffuse log likelihood as follows:

EWMA1fit. <- KFS(EWMA1fit$model)
nInf <- length(EWMA1fit.$Finf)
N <- nrow(EWMA1fit.$model$y)
i. <- (nInf+1):N
-0.5*with(EWMA1fit., sum(log(Finf)) + (N-nInf)*log(2*pi) +
            sum(log(F[i.]) + v[i.]^2/F[i.]) )
## [1] 372.9323

*** JOUNI: I’d still like to manually compute the marginal log(likelihood), and I’m still struggling to figure this out.

We’d also like to compute the one-step ahead root mean square prediction error (rmse). These prediction errors are available as component “v” of the ouput from KFS:

EWMA1.KFS <- KFS(EWMA1fit$model)
sqrt(mean(EWMA1.KFS$v^2))
## [1] 0.04568241

As noted above, log(1+x) is approximately x when x is small. Therefore, we can interpret this as saying that sqrt(mean(EWMA1.v^2)) = 0.046 indicates a root mean square prediction error of approximately 4.6 percent.

Let’s wrap these three goodness-of-fit measures in a fuction, so we can easily compute 3-number summary for subsequent model comparison:

logLik2 <- function(object, ...){
  ll2 <- rep(NA, 4)
  names(ll2) <- c('pars', 'diffuse', 'marginal', 'rmse')
  ll2[1] <- length(object$optim.out$par)
  ll2[2] <- logLik(object$model)
  ll2[3] <- logLik(object$model, marginal=TRUE)
  K. <- KFS(object$model)
  ll2[4] <- sqrt(mean(K.$v^2))
#
  ll2
}

Let’s create a matrix to store these numbers for the different models we’ll fit:

mdls <- c('EWMA1', 'EWMA2', 'EWMA2pres1', 'EWMA2pres2', 
          'EWMA2pres3', 'EWMA2war0', 'EWMA2war1')
logLiks <- matrix(NA, 7, 4, 
  dimnames=list(mdls, 
        c('pars', 'diffuse', 'marginal', 'rmse')) )

Let’s look at and store this “logLik2” summary for this first model:

logLiks[1, ] <- logLik2(EWMA1fit)
logLiks[1, ]
##         pars      diffuse     marginal         rmse 
##   2.00000000 372.93226721 375.64031741   0.04568241

Beyond this, it is also interesting to look at the estimated migration and observation standard deviations:

exp(EWMA1fit$optim.out$par/2)
##          lnQ          lnH 
## 4.578484e-02 2.007486e-08

These estimates put the observation error at essentially zero. It would be interesting to construct a confidence region for these two parameters, but we shall not attempt that here.

To help us understand this model, we compute the one-step ahead predictions with error bounds as follows:

EWMA1pred <- predict(EWMA1fit$model, filtered = TRUE)

Let’s plot this vs. Presidents and wars:

plotPresWars(EWMA1.KFS$v, 
    yHoover=0.8, yFDR=.3, yWar=seq(.6, by=-.1, len=9), 
    ylab='one-step ahead prediction errors')
abline(h=0)
EWMA residuals, Presidents and wars

EWMA residuals, Presidents and wars

There are three negative spikes between 1900 and 1950. Let’s look closer at those:

plotPresWars(EWMA1.KFS$v, xlim=c(1900, 1950), 
    yHoover=0.8, yFDR=.3, yWar=seq(.6, by=-.1, len=9), 
    ylab='one-step ahead prediction errors')
abline(h=0)
EWMA residuals, Presidents and wars

EWMA residuals, Presidents and wars

We understand the Hoover and post-WW II spikes. Let’s take a closer look at the 1908 spike:

plotPresWars(EWMA1.KFS$v, xlim=c(1900, 1910), 
    yHoover=0.8, yFDR=.3, yWar=seq(.6, by=-.1, len=9), 
    ylab='one-step ahead prediction errors')
abline(h=0)
EWMA residuals, Presidents and wars

EWMA residuals, Presidents and wars

Let’s summarize what we get from these plots:

  1. During the Hoover administration, the economy performed consistently worse than this simple Bayesian EWMA predicted. This is what we expected from earlier macroeconomic research, but it’s nice to see it confirmed. (It increases our confidence in the KFAS software.)

  2. During the FDR years including WW II, the economy outperformed the EWMA predictions, except for the Recession of 1937-38. (Again, this increases our confidence in the software.)

  3. We do NOT generally see that other wars are accompanied with large positive residuals, as we would expect if WW II and not FDR ended the Great Depression. There were interesting spikes during the Northwest Indian War, the U.S. Civil War, and the Spanish-American-Philippine war. However, they were isolated and not obviously greater than other spikes that did not accompany wars.

As one more check on this model, let’s set the aberrant period 1929-1947 to NA and refit the model. Let’s first identify this most discrepant period:

depr1 <- with(GDP, (1927<=Year) & (Year <= 1948))
GDP[depr1, c('Year', 'realGDPperCapita')]
##     Year realGDPperCapita
## 171 1927          8267.41
## 172 1928          8259.88
## 173 1929          8669.00
## 174 1930          7847.00
## 175 1931          7288.00
## 176 1932          6308.00
## 177 1933          6192.00
## 178 1934          6817.00
## 179 1935          7373.00
## 180 1936          8273.00
## 181 1937          8643.00
## 182 1938          8292.00
## 183 1939          8881.00
## 184 1940          9583.00
## 185 1941         11171.00
## 186 1942         13138.00
## 187 1943         15165.00
## 188 1944         16181.00
## 189 1945         15850.00
## 190 1946         13869.00
## 191 1947         13456.00
## 192 1948         13776.00

In spite of the stock market crash that began on “black Thursday,” Oct. 24, 1929, the average annual income in 1929 was almost 5 percent higher than in 1928. Average annual income fell in 1930, so we’ll start the exclusion period there. We’ll end the exclusion period in 1947, because the post-war recession ended then, as indicated by the fact that annual income in that year was less than in both 1946 and 1948:

depr <- with(GDP, (1929<=Year)&(Year<=1947))

We can remove these years from any appropriate copy of realGDPperCapita and fit to that modified copy. Conveniently, EWMA1mdl has such a copy in its first component, y:

names(EWMA1mdl)
##  [1] "y"            "Z"            "H"            "T"           
##  [5] "R"            "Q"            "a1"           "P1"          
##  [9] "P1inf"        "u"            "distribution" "tol"         
## [13] "call"         "terms"
EWMA1mdl.woDepr <- EWMA1mdl
EWMA1mdl.woDepr$y[depr] <- NA
EWMA1.woDeprFit<-fitSSM(EWMA1mdl.woDepr, 
                         inits=c(lnQ=0, lnH=0))
EWMA1.woDeprPred <- predict(EWMA1.woDeprFit$model, 
              interval='prediction', filtered=TRUE)

Now let’s plot the predictions:

plotPresWars(GDP$realGDPperCapita/1000, log='y', 
             main=paste0("EWMA fit\n", 
                         'without the Great Depression'), 
             ylab='thousands of 2009 dollars', las=1, 
             type='n', xlab='Year', xlim=c(1920, 1960), 
             ylim=c(5, 20), cex.main=1.1)
text(GDP$realGDPperCapita/1000, 'x', cex=0.5)
matlines(GDP$Year, exp(EWMA1.woDeprPred)/1000, 
         col=c('blue', 'darkgreen', 'darkgreen'), 
         lty=c(1, 2, 2))
EWMA with 1929-1947 set to NA

EWMA with 1929-1947 set to NA

The prediction is flat during the entire period set to NA, as we should expect from the model. Moreover, 1948 is outside the error bounds. This illustrates the need for something more, like a double EWMA. This is a similar model that estimates both level and slope.

2. Double exponential smoothing: More realistic but still not great for the Great Depression

Double exponential smoothing runs an EWMA on the slope of a line and uses that in adjusting the level. The standard Bayesian modification of that adjusts the weight on the last observation in proportion to the information content of the last observation relative to information accumulated in the previous estimate of the state, as suggested above. This can be written using the KFAS notation above as follows:

\[y_t = Z_t \alpha_t + \epsilon_t\]

\[\alpha_{t+1} = T_t \alpha_t + R_t \eta_t.\]

where \(y_t\) is the observation at time \(t\) on log(readGDPperCapita), as before,

\(\alpha_t\) has two components corresponding to the level and slope,

\(Z_t\) is a row vector of (1, 0), so \(y_t\) is the estimated level plus noise,

\(T_t\) is a 2 x 2 matrix with 0 in the lower left and 1’s elsewhere, so the new level is the previous level plus the slope while the new slope is the previous slope, plus migration noise:

\[T_t = \left[\begin{array} {rrr} 1 & 1 \\ 0 & 1 \end{array}\right] \]

Similar to the simple EWMA considered above, \(H_t\) is the observation variance. However, the state vector \(\alpha_t\) is 2-dimensional. To manage this, we assume that \(R_t\) is the 2-dimensional identity matrix, and \(Q_t\) is a 2-dimensional diagonal matrix. This gives us 3 parameters to estimate: The two diagonal elements of the migration variance matrix \(Q_t\) and the observation error variance \(H_t\). As before, by default, the algorithm estimates the logarithms of these parameters, which we will call lnQ.lvl, lnQ.slope, and lnH.

Modeling is a straightforward generalization of the above following the same process:

formula -> SSModel -> fitSSM -> KFS

We start by supplying a formula using SSMtrend(2, …). We select a starting value for \(a_1\) with 0 slope, because that will make the EWMA1 discussed above a special case of this model (with zero migration variance for the slope). That’s important, because the Neyman–Pearson lemma says that the most powerful test comparing two hypotheses is the likelihood ratio, and it helps to have one hypothesis nested within the other, e.g., when the “null hypothesis” is a special case of the “alternative”. Moreover, with most nested hypotheses Wilks’s theorem says that twice the log(likelihood ratio) of nested hypotheses is approximately distributed as chi-square with degrees of freedom equal to number of parameters in the larger model fixed to create the smaller model. Unfortunately, Wilks’s theorem does not apply when the null hypothesis is on a boundary, as with a variance parameter at 0. In that case, the likelihood ratio is still most powerful, but twice the log(likelihood ratio) may not be approximately chi-square. Pinheiro and Bates (2000) recommend simulation to estimate p-values, but provide an example that suggests that 0.5 times the p-value for a chi-square(1) provides a reasonable approximation for the actual p-value when testing if a single variance parameter is zero; we use this approximation below. (They provide other examples that show that this simple rule does not generalize and adjusts in the wrong direction when eliminating fixed effects from a mixed model; see their pp. 87-89.)

As before, we start by specifying a formula for the double EWMA:

(a1.2 <- c(level=log(GDP$realGDPperCapita[1]), slope=0))
##    level    slope 
## 7.009698 0.000000
EWMA2fmla <- log(realGDPperCapita)~
  SSMtrend(2,Q=list(lnQ.lvl=matrix(NA), lnQ.slope=matrix(NA)), 
           a1=a1.2)

As above, we convert this to an SSModel:

EWMA2mdl <-SSModel(EWMA2fmla, GDP, H=matrix(NA))

We pass this as with the simple EWMA to fitSSM:

EWMA2init <- c(lnQ.lvl=0, lnQ.slope=0, lnH=0)
EWMA2fit<-fitSSM(EWMA2mdl, inits=EWMA2init)

We pause here compute the 3 goodless of fit statistics discussed above:

logLiks[2, ] <- logLik2(EWMA2fit)
logLiks[1:2, ]
##       pars  diffuse marginal       rmse
## EWMA1    2 372.9323 375.6403 0.04568241
## EWMA2    3 384.6971 394.2869 0.04274703

The RMSE at 0.043 is only slightly smaller than the root mean square prediction error of 0.046 for the standard Bayesian EWMA discussed above.

Let’s continue with twice the log(likelihood ratio):

(EWMAlogLik2.1 <- 2*diff(logLiks[1:2,1]))
## EWMA2 
##     2

Unfortunately, as discussed above, the assumptions for Wilks’s theorem do not apply here. The standard recommendation would be to compute a p-value for this using Monte Carlo. Before going to that trouble, we first try the modification suggested by Pinheiro and Bates (2000), mentioned above:

0.5*pchisq(EWMAlogLik2.1, 1, lower=FALSE)
##     EWMA2 
## 0.0786496

This gives us 6e-7. Even if this is off by two or three orders of magnitude, the effect is still statistically significant – as we would expect from the plots above.

Let’s continue as above by examining the estimated migration and observation error standard deviations:

exp(EWMA2fit$optim.out$par/2)
##      lnQ.lvl    lnQ.slope          lnH 
## 4.258525e-02 4.808260e-07 2.508964e-24

The migration standard deviation for the level is only slightly smaller that for the EWMA1 above. The migration standard deviation for the slope and measurement noise seem much smaller than we might have expected. It would be interesting to compute a 95 percent confidence interval on the standard deviation of the slope migration; we will not attempt that here. (The preferred approach, as previously indicated, would be Monte Carlo.)

Next, we call KFS:

EWMA2.KFS <- KFS(EWMA2fit$model)

To help us understand this model, we plot the smoothed estimates of the state vector over time:

plotPresWars(EWMA2.KFS$alphahat, yHoover=0.38, yFDR=0.55, 
             yWar=seq(.48, by=-.05, len=9), 
             main=paste0('Double EWMA fit\n', 
                        'to USGDP per capita'), 
             cex.main=1.1)
Estimated level and slope of GDP

Estimated level and slope of GDP

The slope has the shape we expected from the plot, but it does not seem to fit the Hoover and FDR years well. Because of the poor labeling of the vertical axis, let’s compute the ranges of the smoothed state estimates:

apply(EWMA2.KFS$alphahat, 2, range)
##          level      slope
## [1,]  7.009698 0.01704459
## [2,] 10.827687 0.01704460

The estimated slope is surprisingly constant, ranging only from 0.01704459 to 0.01704460. That’s much stiffer than might have been expected from other analyses that suggest that the growth in the early nineteenth century was less than in the late twentieth. To check this, we’ll plot 20-year running averages of the growth:

avgGrowth20 <- caTools::runmean(
        diff(log(GDP$realGDPperCapita)), 20)
plotPresWars(GDP$Year[-1], avgGrowth20, type='l', 
             yFDR=0.3, las=1, xlab='Year', 
             ylab='20-year avg annual % growth', 
             main=paste0('Smoothed avg growth\n',
                    'in annual income per capita') )
abline(h=c(0, .02, .04))
20-year running averages of the growth rate

20-year running averages of the growth rate

One might naively expect to see variability in alphahat[, ‘slope’] closer to what we see in this plot. Evidently, the maximum likelihood smoothing window is closer to 200 years than 20. It would be helpful to examine profile likelihood plots over the smallest two estimated parameters, lnQ.slope and lnH, as discussed, e.g., with the help page for plot.profile.nls and in the lme4 vignette (Bates, Mächler, Bolker, and Walke, 2014) and see how the above plot of alphahat[, ‘slope’] would change when increasing lnQ.slope to its upper 80 percent confidence limit. We shall not pursue that here.

As another evaluation of the fit during the FDR years, let’s plot the one-step ahead prediction errors as before:

plotPresWars(EWMA2.KFS$v, type='l', 
     main='Double EWMA\nprediction errors', 
     ylab='log(2009 dollars)', 
     yHoover=.8, yFDR=0.1, 
             yWar=seq(.7, by=-.12, len=9))
Double EWMA, Presidents and wars

Double EWMA, Presidents and wars

Before proceeding with the next model, let’s repeat the analysis with the Great Depression set to NA:

EWMA2mdl.woDepr <- EWMA2mdl
EWMA2mdl.woDepr$y[depr] <- NA
EWMA2.woDeprFit<-fitSSM(EWMA2mdl.woDepr, inits=EWMA2init)
EWMA2.woDeprPred <- predict(EWMA2.woDeprFit$model, 
              interval='prediction', filtered=TRUE)
plotPresWars(GDP$realGDPperCapita/1000, log='y', 
             main=paste0('Double EWMA fit\n', 
                         'without the Great Depression'), 
             ylab='thousands of 2009 dollars', las=1, 
             type='n', xlab='Year', xlim=c(1920, 1960), 
             ylim=c(5, 20), cex.main=1.1) 
text(GDP$realGDPperCapita/1000, 'x', cex=0.5)
matlines(GDP$Year, exp(EWMA2.woDeprPred)/1000, 
         col=c('blue', 'darkgreen', 'darkgreen'), 
         lty=c(1, 2, 2))
matlines(exp(EWMA2.woDeprPred)/1000, 
         col=c('blue', 'darkgreen', 'darkgreen'), 
         lty=c(1, 2, 2))
Double EWMA with 1929-1947 set to NA

Double EWMA with 1929-1947 set to NA

The predictions are better than the standard EWMA, but the prediction limits are still too tight: Most of Hoover’s presidency is below the lower prediction limit. FDR started below the lower limit and ended above the upper limit. Clearly this period cannot be described by the model fit to the rest of the data.

Let’s look at the growth numbers for this period:

sel27.49 <- (GDP$Year %in% 1927:1949)
cbind(1928:1949, diff(EWMA2.woDeprPred[sel27.49,'fit']))
##       [,1]         [,2]
##  [1,] 1928 -0.004384176
##  [2,] 1929 -0.001024186
##  [3,] 1930  0.014561401
##  [4,] 1931  0.014561401
##  [5,] 1932  0.014561401
##  [6,] 1933  0.014561401
##  [7,] 1934  0.014561401
##  [8,] 1935  0.014561401
##  [9,] 1936  0.014561401
## [10,] 1937  0.014561401
## [11,] 1938  0.014561401
## [12,] 1939  0.014561401
## [13,] 1940  0.014561401
## [14,] 1941  0.014561401
## [15,] 1942  0.014561401
## [16,] 1943  0.014561401
## [17,] 1944  0.014561401
## [18,] 1945  0.014561401
## [19,] 1946  0.014561401
## [20,] 1947  0.014561401
## [21,] 1948  0.014561401
## [22,] 1949  0.236245855

Without the data for 1930 to 1948, the model predicts an increase of almost 1.5 percent per year (exp(0.01456)-1 = 1.467 percent) for each year with no data, then a jump of 0.236 in 1949 to close the gap created by the faster-than-average growth during the FDR years, even ignoring the depth of the depression in the year that Hoover left office.

Overall, this is a clear improvement over the standard EWMA discussed above.

Let’s see if we can do better adding a random effect for presidents. We consider five different ways to do this, as outlined in the table above.

3. Allowing the slope to change only between executives

A simple presidential effect model is to assume that the slope is constant within each presidency but can change between administrations. This uses the same state transition equation as the double EWMA just discussed, but the migration variance is zero except for the last year of each presidency: That allows the slope to jump for the first year of the next presidency.

This model can be estimated using KFAS but requires using a state transition variance, Q, that is not constant. To fit this, we must write a special “updatefn” for fitSSM, because the default “updatefn” will not handle a time varying Q:

(make the migration variance, Q, time varying) -> formula -> SSModel

-> (create updatefn) -> fitSSM -> KFS

The state transition variance, Q, is typically specified as a list in SSMtrend. Before changing this, let’s first look at the structure of EWMA2mdl:

str(EWMA2mdl)
## List of 14
##  $ y           : ts [1:225, 1] 7.01 7.04 7.08 7.13 7.22 ...
##   ..- attr(*, "tsp")= num [1:3] 1790 2014 1
##  $ Z           : num [1, 1:2, 1] 1 0
##   ..- attr(*, "dimnames")=List of 3
##   .. ..$ : NULL
##   .. ..$ : chr [1:2] "level" "slope"
##   .. ..$ : NULL
##  $ H           : logi [1, 1, 1] NA
##  $ T           : num [1:2, 1:2, 1] 1 0 1 1
##   ..- attr(*, "dimnames")=List of 3
##   .. ..$ : chr [1:2] "level" "slope"
##   .. ..$ : chr [1:2] "level" "slope"
##   .. ..$ : NULL
##  $ R           : num [1:2, 1:2, 1] 1 0 0 1
##   ..- attr(*, "dimnames")=List of 3
##   .. ..$ : chr [1:2] "level" "slope"
##   .. ..$ : NULL
##   .. ..$ : NULL
##  $ Q           : num [1:2, 1:2, 1] NA 0 0 NA
##  $ a1          : num [1:2, 1] 7.01 0
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:2] "level" "slope"
##   .. ..$ : NULL
##  $ P1          : num [1:2, 1:2] 0 0 0 0
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:2] "level" "slope"
##   .. ..$ : chr [1:2] "level" "slope"
##  $ P1inf       : num [1:2, 1:2] 1 0 0 1
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:2] "level" "slope"
##   .. ..$ : chr [1:2] "level" "slope"
##  $ u           : chr "Omitted"
##  $ distribution: chr "gaussian"
##  $ tol         : num 1.49e-08
##  $ call        : language SSModel(formula = EWMA2fmla, data = GDP, H = matrix(NA))
##  $ terms       :Classes 'terms', 'formula' length 3 log(realGDPperCapita) ~ SSMtrend(2, Q = list(lnQ.lvl = matrix(NA),      lnQ.slope = matrix(NA)), a1 = a1.2)
##   .. ..- attr(*, "variables")= language list(log(realGDPperCapita), SSMtrend(2, Q = list(lnQ.lvl = matrix(NA),      lnQ.slope = matrix(NA)), a1 = a1.2))
##   .. ..- attr(*, "factors")= int [1:2, 1] 0 1
##   .. .. ..- attr(*, "dimnames")=List of 2
##   .. .. .. ..$ : chr [1:2] "log(realGDPperCapita)" "SSMtrend(2, Q = list(lnQ.lvl = matrix(NA), lnQ.slope = matrix(NA)), a1 = a1.2)"
##   .. .. .. ..$ : chr "SSMtrend(2, Q = list(lnQ.lvl = matrix(NA), lnQ.slope = matrix(NA)), a1 = a1.2)"
##   .. ..- attr(*, "term.labels")= chr "SSMtrend(2, Q = list(lnQ.lvl = matrix(NA), lnQ.slope = matrix(NA)), a1 = a1.2)"
##   .. ..- attr(*, "specials")=Dotted pair list of 6
##   .. .. ..$ SSMregression: NULL
##   .. .. ..$ SSMtrend     : int 2
##   .. .. ..$ SSMseasonal  : NULL
##   .. .. ..$ SSMcycle     : NULL
##   .. .. ..$ SSMarima     : NULL
##   .. .. ..$ SSMcustom    : NULL
##   .. ..- attr(*, "order")= int 1
##   .. ..- attr(*, "intercept")= int 1
##   .. ..- attr(*, "response")= int 1
##   .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
##  - attr(*, "class")= chr "SSModel"
##  - attr(*, "p")= int 1
##  - attr(*, "m")= int 2
##  - attr(*, "k")= int 2
##  - attr(*, "n")= int 225
##  - attr(*, "tv")= int [1:5] 0 0 0 0 0
##  - attr(*, "state_types")= chr [1:2] "level" "slope"
##  - attr(*, "eta_types")= chr [1:2] "level" "slope"

The most important thing to note for present purposes is that the component Q is “num [1:2, 1:2, 1] NA 0 0 NA”:

EWMA2mdl$Q
## , , 1
## 
##      [,1] [,2]
## [1,]   NA    0
## [2,]    0   NA

This was specified above as “Q=list(lnQ.lvl=matrix(NA), lnQ.slope=matrix(NA))” in the call above to SSMtrend that created EWMA2fmla. We will change that here by replacing “lnQ.slope=matrix(NA)” with an array of all zeros except for making the slope component NA for the last year of each administration:

EWMA2pres1Q.slope <- array(0, c(1,1,N), 
      dimnames=list('slope', 'slope', GDP$Year))
contin <- duplicated(GDP$executive)
presLastYr <- c(!contin[-1], FALSE)
EWMA2pres1Q.slope[,,presLastYr] <- NA

We now create a formula for this model:

EWMA2pres1fmla <- log(realGDPperCapita) ~ 
  SSMtrend(2, Q = list(lnQ.lvl=matrix(NA), 
                       lnQ.slope=EWMA2pres1Q.slope), a1=a1.2)

As before, we use this to build an SSModel:

EWMA2pres1mdl <-SSModel(EWMA2pres1fmla, GDP, H=matrix(NA))

The next step is to create the desired updatefn assuming we estimate three log(variance) parameters for lnQ.lvl, lnQ.slope, and lnH, as for the double EWMA above:

EWMA2pres1updt <- function(pars=EWMA2init, model=EWMA2pres1mdl, 
                           ...){
  vars <- exp(pars)
  Q <- model$Q
  Q[1,1,] <- vars[1]
  Q[2,2,presLastYr] <- vars[2]
  model$Q <- Q
  model$H[] <- vars[3]
  model
}

We’re now ready to call fitSSM:

EWMA2pres1fit <- fitSSM(EWMA2pres1mdl, inits=EWMA2init, 
                        updatefn=EWMA2pres1updt)
logLiks[3, ] <- logLik2(EWMA2pres1fit)
logLiks[1:3, ]
##            pars  diffuse marginal       rmse
## EWMA1         2 372.9323 375.6403 0.04568241
## EWMA2         3 384.6971 394.2869 0.04274703
## EWMA2pres1    3 384.6971 394.2868 0.04274704

Both log(likelihoods) and the rmse are virtually identical to those numbers for the double EWMA above.

As before, we’ll compute and plot the smoothed estimate of the state vector:

EWMA2pres1.KFS <- KFS(EWMA2pres1fit$model)
plotPresWars(EWMA2pres1.KFS$alphahat, yHoover=0.38, yFDR=0.59, 
              yWar=seq(.48, by=-.05, len=9), 
              main=paste0('Smoothed level and slope\n', 
                          'for first President model'), 
             cex.main=1.1)
Estimated level and slope of GDP for first President model

Estimated level and slope of GDP for first President model

This is very similar to the double EWMA but with the slope held constant within administration – precisely what we attempted to fit.

However, since the model seems not to improve upon the double EWMA, we will leave this here and try the second model mentioned above.

4. Assuming no continuity in slope between executives

Another approach to modeling a president effect is to restart the slope from 0 at the beginning of each presidency. We will do this by modifying the EWMA2pres1 model above so the slope component of the state transition matrix is 0 for the last year of each presidency and 1 otherwise:

\[T_t = \left[\begin{array} {rrr} 1 & 1 \\ 0 & \gamma_t \end{array}\right], \]

where \(\gamma_t\) = 0 for the last year of a presidency, and 1 otherwise. This makes \(\alpha_{2,t+1} = \eta_t\) for \(t\) = the last year of each presidency (and \(t+1\) = the first year of the next). Then with \(Q_{2,2,t}\) = 0 except for the last year of each presidency, that means that \(\alpha_{2,t+1} = \alpha_{2,t}\) for continuing years.

We’ll accomplish this by adding another step to the process used above, modifying the SSModel before fitting it:

(make the migration variance, Q, time varying) -> formula -> SSModel ->

(make the state transition matrix, T, time varying) -> (create updatefn) -> fitSSM -> KFS

In this particular case, we can use the same Q and formula as in the previous section, and modify the previous model to make T time varying as desired:

EWMA2pres2mdl <- EWMA2pres1mdl
T.other <- EWMA2pres1mdl$T
EWMA2pres2.T <- array(T.other, c(2,2,N), 
      dimnames=dimnames(T.other))
dimnames(EWMA2pres2.T)[[3]] <- GDP$Year
EWMA2pres2.T[2,2,presLastYr] <- 0
EWMA2pres2mdl$T <- EWMA2pres2.T

We’re not quite done modifying this model: As noted with help(‘SSModel’), we must set attribute ‘tv’ to reflect that both T and Q are time varying. From “str(EWMA2mdl)” above, we see that attr(EWMA2mdl1, “tv”)= int [1:5] 0 0 0 0 0, indicating a model in which none of Z, H, T, R, and Q are time varying. This is different for EWMA2pres1mdl:

attr(EWMA2pres1mdl, 'tv')
## [1] 0 0 0 0 1

This says that only Q is time varying. We just made T time varying and must change this ‘tv’ attribute to make it official. In the process, we will add names to make the structure easier to read:

attr(EWMA2pres2mdl, 'tv') <- c(Z=0L, H=0L, T=1L, R=0L, Q=1L)

Formally at least, we have the same parameters to estimate as with EWMA2pres1. However, we can set \(P_{2,2,1}\) to match \(Q_{2,2,t}\) for the last year of each presidency; note that P1 = “num [1:2, 1:2, 1] 0 0 0 0” in “str(EWMA2mdl)” above.

EWMA2pres2updt <- function(pars=EWMA2init, model=EWMA2pres1mdl, 
                           ...){
  vars <- exp(pars)
  Q <- model$Q
  Q[1,1,] <- vars[1]
  Q[2,2,presLastYr] <- vars[2]
  model$Q <- Q
  diag(model$P1) <- vars[1:2]
  model$H[] <- vars[3]
  model
}

We now fit this model as follows:

{r} EWMA2pres2fit <- fitSSM(EWMA2pres2mdl, inits=EWMA2init, updatefn=EWMA2pres2updt) logLiks[4, ] <- logLik2(EWMA2pres2fit) logLiks[1:4, ]

This has a slightly better fit, as indicated by the increase in the log(likelihoods) and reduction in rmse.

We next compute and plot the smoothed estimate of the state vector, as before:

{r,fig.cap = "Estimated level and slope of GDP, decoupling the presidents"} EWMA2pres2.KFS <- KFS(EWMA2pres2fit$model) plotPresWars(EWMA2pres2.KFS$alphahat, yHoover=0.72, yFDR=.3, yWar=seq(.6, by=-.08, len=9), main='Decoupled presidents')

The slope here shows substantially more variability. Let’s look at it more closely:

{r,fig.cap = "Estimated level slope of GDP, decoupling the presidents"} plotPresWars(EWMA2pres2.KFS$alphahat[, 'slope'], yHoover=0.72, yFDR=.3, yWar=seq(.6, by=-.08, len=9), ylab='slope') abline(h=0)

This strongly suggests that the Hoover administration was a statistical outlier, while FDR (with World War II) may or may not have been. To investigate this further, let’s look at the slope and it’s variance for each president. help(‘KFS’) tells us that the “V” component is the covariance of “alphahat”, given all the data (and the “P” and “Pinf” components are the non-diffuse and diffuse parts of the error covariance of the one-step ahead predicted states in “a”). Just to make sure we understand this correctly, let’s plot “V”:

{r,fig.cap = "Variance of the estimated slope, decoupling the presidents"} plotPresWars(GDP$Year, EWMA2pres2.KFS$V[2,2,], yHoover=0.72, yFDR=.3, yWar=seq(.6, by=-.08, len=9), type='l', ylab='smoothed slope variance')

This is plausibly constant within administration, which we would expect from specifying zero migration variance for the slope. If this is correct, we would expect the “information” (i.e., reciprocal variance) to be proportional to the number of years in each presidency. To confirm that, we compute the posterior variance and the number of years for each presidency:

{r,fig.cap = "Slope, variance, and number of years for each presidency"} EWMA2pres2V <- tapply(EWMA2pres2.KFS$V[2,2,], GDP$executive, mean) EWMA2pres2n <- table(GDP$executive) plot(as.numeric(EWMA2pres2n), 1/EWMA2pres2V, xlab='Number of years', ylab='info = 1/(smoothed variance)')

This looks like a straight line except for two points. Who are those two?

{r,fig.cap = "Posterior variance vs. number of years for each presidency, labeled"} plot(as.numeric(EWMA2pres2n), 1/EWMA2pres2V, xlab='Number of years', ylab='info = 1/(smoothed variance)') text(as.numeric(EWMA2pres2n), 1/EWMA2pres2V, 1:41)

Of course: The two off the line are the first and last.

Let’s compute the slope and 95 percent confidence limits:

{r,fig.cap = "Decoupled slope for each presidency"} EWMA2pres2pres <- tapply(EWMA2pres2.KFS$alphahat[, 'slope'], GDP$executive, mean) sz2 <- qnorm(.975)*sqrt(EWMA2pres2V) sz2. <- sz2 %o% c(lower=-1, upper=1) EWMA2pres2slim <- as.numeric(EWMA2pres2pres) + sz2. EWMA2pres2Slim <- cbind(data.frame( executive=ordered(names(sz2), names(sz2)), slope=EWMA2pres2pres), EWMA2pres2slim) library(ggplot2) ggplot(EWMA2pres2Slim, aes(x=executive, y=slope))+geom_point() + coord_cartesian(ylim = c(-0.08, 0.08)) + geom_errorbar(ymin=EWMA2pres2Slim$lower, ymax=EWMA2pres2Slim$upper)+ geom_hline(aes(yintercept=mean(slope),color="red")) # mean slope

This clearly says that the Hoover and FDR administrations were different from all others, but it still doesn’t say whether FDR himself was different or World War II created the difference. Before considering that, let’s continue with the other alternative presidency model.

5. EWMA2 with an adjustment for each president

This model considers a 3-dimensional state vector: level, slope and a Presidential adjustment for slope. For this, we will want a state transition matrix like the following:

\[T_t = \left[\begin{array} {rrr} 1 & 1 & 1 \\ 0 & 1 & 0 \\ 0 & 0 & \gamma_t \end{array}\right], \]

where \(\gamma_t\) = 0 for the last year of a presidency, and 1 otherwise, as with immediately preceding model, EWMA2pres2.

We will create this by modifying a triple EWMA following a process similar to the above:

(make the migration variance, Q, time varying) -> formula -> SSModel ->

(make the state transition matrix, T, time varying) -> (create updatefn) -> fitSSM -> KFS

Let’s start by modifying the formula, similar to what we did with EWMA2pres1fmla, above, but for the new 3-dimensional state

a1.3 <- c(level=log(GDP$realGDPperCapita[1]), slope=0, 
           Pres=0)
EWMA2pres3fmla <- log(realGDPperCapita) ~ 
  SSMtrend(3, Q = list(level=matrix(NA), slope=matrix(NA),
                       Pres=EWMA2pres1Q.slope),
           a1=a1.3)

We next create an SSModel, as before:

EWMA2pres3mdl <- SSModel(EWMA2pres3fmla, GDP, H=matrix(NA))

Now we need to create the desired time-varying state transition matrix, T:

pres3states <- names(a1.3)
pres3T <- array(EWMA2pres3mdl$T, dim=c(3,3,N), 
        dimnames=list(pres3states, pres3states, GDP$Year) )
pres3T[1,3,] <- 1
pres3T[2,3,] <- 0 
pres3T[3,3, presLastYr] <- 0
EWMA2pres3mdl$T <- pres3T

As for EWMA2pres2mdl, we need also to fix the ‘tv’ attribute:

attr(EWMA2pres3mdl, 'tv') <- c(Z=0L, H=0L, T=1L, R=0L, Q=1L)

We now have 4 parameters to estimate:

EWMA2pres3init <- c(lnQ.lvl=0, lnQ.slope=0, lnQ.Pres=0, lnH=0)

And the state vector now has 3 dimensions, not 2 as before. Both these changes require a different ‘updatefn’:

EWMA2pres3updt <- function(pars=EWMA2pres3init, 
                      model=EWMA2pres3mdl, ...){
  vars <- exp(pars)
  Q <- model$Q
  Q[1,1,] <- vars[1]
  Q[2,2,] <- vars[2]
  Q[3,3,presLastYr] <- vars[3]
  model$Q <- Q
  model$H[] <- vars[4]
  model
} 
EWMA2pres3fit <- fitSSM(EWMA2pres3mdl, inits=EWMA2pres3init, 
                        updatefn=EWMA2pres3updt) 
logLiks[5, ] <- logLik2(EWMA2pres3fit)
logLiks[1:5, ]
##            pars  diffuse marginal       rmse
## EWMA1         2 372.9323 375.6403 0.04568241
## EWMA2         3 384.6971 394.2869 0.04274703
## EWMA2pres1    3 384.6971 394.2868 0.04274704
## EWMA2pres2   NA       NA       NA         NA
## EWMA2pres3    4 391.4337 403.4437 0.04106114

This is the best logLik and rmse so far. However, the increase in logLik is slightly less than 2, which makes the AIC and the BIC worse.

NOTE: An equivalent model estimating the same parameters could have a 43-dimensional state vector with components level, slope, and Pres1, …, Pres41. The memory and compute time required to fit it would be greater, but the logLik and parameter estimates would be the same. This is unnecessary here but could be required with, e.g., an experiment with levels for a random effect that are not sequential in time.

6. EWMA2 with a fixed effect for war

The plots discussed above suggest that if war has a macroeconomic effect, it increases the rate of economic growth. This suggests we model the impact of war as an adjustment to slope for incrementing level within the state transition matrix, \(T_t\). It further suggests we want a four-dimensional state vector, with one component (\(w_0\)) for the presence (\(w_0\) = 1) or absence (\(w_0\) = 0) of war and the second proportional to battle deaths per million population (\(w_1\)):

\[T_t = \left[\begin{array} {rrrr} 1 & 1 & w_0 & w_1 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{array}\right], \]

For this model, Q is constant. Therefore, we do not need to create a time-varying Q nor an updatefn as we did for models 3-5. This simplifies the modeling process to the following:

formula -> SSModel ->

(make the state transition matrix, T, time varying) -> fitSSM -> KFS

We can model war as a fixed effect multiple ways. One is to use \(R_t\) = the 4-dimensional identity matrix as before, with \(Q_t\) being a diagonal matrix with only the first two elements nonzero. Another is to let \(R_t\) be the first two columns of a 4-dimensional identity matrix and estimate the diagonal of the resulting 2 x 2 \(Q_t\) matrix. We’ll use the former here, because the next model with war as a random effect requires \(R_t\) = the 4-dimensional identity matrix.

We’ll use a similar formula to those above, requesting a fourth degree SSMtrend and modifying \(T_t\) as before:

a1.4 <- c(level=log(GDP$realGDPperCapita[1]), slope=0, 
           war0=0, war1=0)
EWMA2war0fmla <- log(realGDPperCapita) ~ 
    SSMtrend(4, Q = list(level=matrix(NA), slope=matrix(NA),
                       war0=matrix(0), war1=matrix(0)),
           a1=a1.4) 
EWMA2war0mdl <- SSModel(EWMA2war0fmla, GDP, H=matrix(NA))

Now we need to create the desired time-varying state transition matrix, \(T_t\):

warStates <- names(a1.4)
war0T <- array(diag(4), dim=c(4,4,N), 
    dimnames=list(warStates, warStates, GDP$Year))
war0T[1, 2,] <- 1  
war0T[1, 3,] <- (GDP$war!='')
war0T[1, 4,] <- GDP$battleDeathsPMP

Let’s check this \(T_t\) with years at war and peace.

GDP[c(1, 2, 7), c('Year', 'war', 'battleDeathsPMP')]
##    Year                  war battleDeathsPMP
## 34 1790 Northwest Indian War           29.49
## 35 1791 Northwest Indian War           28.62
## 40 1796                                 0.00
war0T[,,c(1, 2, 7)]
## , , 1790
## 
##       level slope war0  war1
## level     1     1    1 29.49
## slope     0     1    0  0.00
## war0      0     0    1  0.00
## war1      0     0    0  1.00
## 
## , , 1791
## 
##       level slope war0  war1
## level     1     1    1 28.62
## slope     0     1    0  0.00
## war0      0     0    1  0.00
## war1      0     0    0  1.00
## 
## , , 1796
## 
##       level slope war0 war1
## level     1     1    0    0
## slope     0     1    0    0
## war0      0     0    1    0
## war1      0     0    0    1

They match. Let’s save this T as part of EWMA2war0mdl:

EWMA2war0mdl$T <- war0T

We also need to fix the ‘tv’ attribute, as we did for models 4 and 5 (EWMA2pres2 and EWMA2pres3):

attr(EWMA2war0mdl, 'tv') <- c(Z=0L, H=0L, T=1L, R=0L, Q=0L)

Let’s also change the names of EWMA2war0mdl$a1 from the default to warStates:

dimnames(EWMA2war0mdl$a1)[[1]] <- warStates

This model has the same 3 parameters to estimate as models 2, 3 and 4 above: lnQ.lvl, lnQ.slope, and lnH. We will use the same initial values: EWMA2init.

(EWMA2war0fit <- fitSSM(EWMA2war0mdl, inits=EWMA2init)) 
## $optim.out
## $optim.out$par
##    lnQ.lvl  lnQ.slope        lnH 
##  -6.312221 -25.341160 -96.748732 
## 
## $optim.out$value
## [1] -370.9859
## 
## $optim.out$counts
## function gradient 
##      154       NA 
## 
## $optim.out$convergence
## [1] 0
## 
## $optim.out$message
## NULL
## 
## 
## $model
## Call:
## SSModel(formula = EWMA2war0fmla, data = GDP, H = matrix(NA))
## 
## State space model object of class SSModel
## 
## Dimensions:
## [1] Number of time points: 225
## [1] Number of time series: 1
## [1] Number of disturbances: 4
## [1] Number of states: 4
## Names of the states:
## [1]  level  slope  war0   war1 
## Distributions of the time series:
## [1]  gaussian
## 
## Object is a valid object of class SSModel.
logLiks[6, ] <- logLik2(EWMA2war0fit)
logLiks[1:6, ]
##            pars  diffuse marginal       rmse
## EWMA1         2 372.9323 375.6403 0.04568241
## EWMA2         3 384.6971 394.2869 0.04274703
## EWMA2pres1    3 384.6971 394.2868 0.04274704
## EWMA2pres2   NA       NA       NA         NA
## EWMA2pres3    4 391.4337 403.4437 0.04106114
## EWMA2war0     3 370.9859 394.7362 0.04644578

Conclusion:

  1. If we use the marginal likelihood rather than the diffuse likelihood, it says that EWMA2war0 gives a slightly better fit: 394.7362 vs. 394.2869. In other words, the difference was primarily due to the use of diffuse likelihood, which is different with the addition of another component to the state vector.

  2. However, the improvement in the fit is so slight, it’s clearly not statistically significant.

Let’s try setting the war components of P1inf to 0 to see if that forces the model to match EWMA2:

EWMA2war0mdl2 <- EWMA2war0mdl 
war0P1inf <- diag(c(1, 1, 0, 0))
dimnames(war0P1inf) <- list(warStates, warStates)
EWMA2war0mdl2$P1inf <- war0P1inf
EWMA2war0fit2 <- fitSSM(EWMA2war0mdl2, inits=EWMA2init)
logLik2(EWMA2war0fit2)
##         pars      diffuse     marginal         rmse 
##   3.00000000 384.69714345 394.28688105   0.04274703

This gives the same logLik and rmse as EWMA2, as we might expect. Let’s modify the modeling process to include an updatefn to allow us to estimate P1[3:4, 3:4] different from diag(0, 2):

formula -> SSModel -> (make the state transition matrix, T, time varying)

(create updatefn) -> fitSSM -> KFS

P1[3:4, 3:4] is a 2 x 2 positive definite matrix. A standard representation for such a matrix is provided by nlme::pdLogChol. The pdLogChol representation of a symmetric positive definite matrix is a vector starting with the logarithms of the diagonal of the Choleski factorization of that matrix followed by its upper triangular portion. To see this, consider the following simple example:

(pd4 <- nlme::pdLogChol(1:6))
## Positive definite matrix structure of class pdLogChol representing
##           [,1]     [,2]      [,3]
## [1,]  7.389056 10.87313  13.59141
## [2,] 10.873127 70.59815  64.33434
## [3,] 13.591409 64.33434 464.42879
(pd4c <- chol(pd4))
##          [,1]     [,2]     [,3]
## [1,] 2.718282 4.000000  5.00000
## [2,] 0.000000 7.389056  6.00000
## [3,] 0.000000 0.000000 20.08554
log(diag(pd4c))
## [1] 1 2 3
pd4c[upper.tri(pd4c)]
## [1] 4 5 6

We use this is the updatefn:

EWMA2war0init <- c(EWMA2fit$optim.out$par, 
                lnP1.33=0, lnP1.44=0, lnP1.34=0)
EWMA2war0updt <- function(pars=EWMA2war0init, 
                      model=EWMA2war0mdl2, ...){
  vars <- exp(pars[1:3])
  Q <- model$Q
  Q[1,1,] <- vars[1]
  Q[2,2,] <- vars[2]
  model$Q <- Q
  model$H[] <- vars[3]
  P1 <- nlme::pdLogChol(pars[4:6])
  model$P1[3:4, 3:4] <- as.matrix(P1)
  model
} 
EWMA2war0fit3 <- fitSSM(EWMA2war0mdl2, inits=EWMA2war0init, 
                        updatefn=EWMA2war0updt)
logLik2(EWMA2war0fit3)
##         pars      diffuse     marginal         rmse 
##   6.00000000 382.62082072 392.21055832   0.04294742

This is worse than EWMA2war0fit and EWMA2war0fit2. Let’s check with starting values designed to give the same answer as EWMA2:

EWMA2war0init2 <- c(EWMA2fit$optim.out$par, 
                lnP1.33=-999, lnP1.44=-999, lnP1.34=0)
EWMA2war0fit4 <- fitSSM(EWMA2war0mdl2, inits=EWMA2war0init2, 
                        updatefn=EWMA2war0updt)
logLik2(EWMA2war0fit4)
##         pars      diffuse     marginal         rmse 
##   6.00000000 384.69714345 394.28688105   0.04274703

This matches logLik(EWMA2war0fit2$model), as expected. Let’s try something between EWMA2war0init and EWMA2war0init2:

EWMA2war0init3 <- c(EWMA2fit$optim.out$par, 
                lnP1.33=-9, lnP1.44=-9, lnP1.34=0)
EWMA2war0fit5 <- fitSSM(EWMA2war0mdl2, inits=EWMA2war0init3,
                        updatefn=EWMA2war0updt)
logLik2(EWMA2war0fit5)
##         pars      diffuse     marginal         rmse 
##   6.00000000 383.82393477 393.41367237   0.04282428

At logLik = 383.8239, this is not as good as EWMA2war0fit5. We might get something better with different starting values.

EWMA2war0init4 <- c(EWMA2fit$optim.out$par, 
                lnP1.33=-2, lnP1.44=-9, lnP1.34=0)
EWMA2war0fit6 <- fitSSM(EWMA2war0mdl2, inits=EWMA2war0init4, 
                        updatefn=EWMA2war0updt)
logLik2(EWMA2war0fit6)
##         pars      diffuse     marginal         rmse 
##   6.00000000 381.54416578 391.13390338   0.04300203

Let’s give up for the moment.

7. EWMA2 with a random effect for war

This model combines EWMA2pres3 with EWMA2war0. We will use a time-varying state transition matrix as follows:

\[T_t = \left[\begin{array} {rrrr} 1 & 1 & w_0 & w_1 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & \omega_t & 0 \\ 0 & 0 & 0 & \omega_t \end{array}\right], \]

where \(\omega_t\) = 1 during wars and 0 otherwise. And we need to make the migration variance time varying to match:

(make the migration variance, Q, time varying) -> formula -> SSModel ->

(make the state transition matrix, T, time varying) -> (create updatefn) -> fitSSM -> KFS

Similar to EWMA2pres2 and EWMA2pres3, we will allow Q for the (war0, war1) components of the state vector to be nonzero the year before the start of each war:

(beforeWar <- which(c(
      head(GDP$war, -1)=='' & tail(GDP$war, -1)!=''), FALSE))
## [1]  22  56  71 108 127 151 160 175

Let’s check:

GDP[beforeWar[1]+0:1,]
##    Year   CPI GDPdeflator population.K realGDPperCapita executive
## 55 1811 12.73        6.39         7436          1614.19   Madison
## 56 1812 12.89        6.30         7651          1631.42   Madison
##            war battleDeaths battleDeathsPMP Keynes
## 55                     0.00            0.00      0
## 56 War of 1812       444.06           58.04      0

We want Q[1:2, 1:2, ] to be diag(NA, 2), Q[3:4, 3:4, beforeWar] to NA; all other elements of Q should be 0:

EWMA2warQ <- array(0, c(4,4,N), 
    dimnames=list(warStates, warStates, GDP$Year))
EWMA2warQ[1,1,] <- NA
EWMA2warQ[2,2,] <- NA
EWMA2warQ[3:4,3:4,beforeWar] <- NA

Let’s confirm we got what we intended:

EWMA2warQ[,,beforeWar[1]+(-1):1]
## , , 1810
## 
##       level slope war0 war1
## level    NA     0    0    0
## slope     0    NA    0    0
## war0      0     0    0    0
## war1      0     0    0    0
## 
## , , 1811
## 
##       level slope war0 war1
## level    NA     0    0    0
## slope     0    NA    0    0
## war0      0     0   NA   NA
## war1      0     0   NA   NA
## 
## , , 1812
## 
##       level slope war0 war1
## level    NA     0    0    0
## slope     0    NA    0    0
## war0      0     0    0    0
## war1      0     0    0    0

SSMtrend(4, …) requires Q to be a list of length 4, so we can’t use that. Instead, we’ll copy the previous SSModel and modify Q directly:

EWMA2war1mdl <- EWMA2war0mdl
EWMA2war1mdl$Q <- EWMA2warQ
attr(EWMA2war1mdl, 'tv')[5] <- 1L

We also want to set T[3:4, 3:4, GDP$war==’’] to 0:

war1T <- war0T
war1T[3:4, 3:4, GDP$war==''] <- 0
war1T[,,c(1:2, 8, 22:24)]
## , , 1790
## 
##       level slope war0  war1
## level     1     1    1 29.49
## slope     0     1    0  0.00
## war0      0     0    1  0.00
## war1      0     0    0  1.00
## 
## , , 1791
## 
##       level slope war0  war1
## level     1     1    1 28.62
## slope     0     1    0  0.00
## war0      0     0    1  0.00
## war1      0     0    0  1.00
## 
## , , 1797
## 
##       level slope war0 war1
## level     1     1    0    0
## slope     0     1    0    0
## war0      0     0    0    0
## war1      0     0    0    0
## 
## , , 1811
## 
##       level slope war0 war1
## level     1     1    0    0
## slope     0     1    0    0
## war0      0     0    0    0
## war1      0     0    0    0
## 
## , , 1812
## 
##       level slope war0  war1
## level     1     1    1 58.04
## slope     0     1    0  0.00
## war0      0     0    1  0.00
## war1      0     0    0  1.00
## 
## , , 1813
## 
##       level slope war0   war1
## level     1     1    1 104.58
## slope     0     1    0   0.00
## war0      0     0    1   0.00
## war1      0     0    0   1.00

Looks good. Let’s put it in the model:

EWMA2war1mdl$T <- war1T

We also need a matching updatefn:

{r} EWMA2war1init <- c(EWMA2fit$optim.out$par, lnP1.33=0, lnP1.44=0, lnP1.34=0) EWMA2war1updt <- function(pars=EWMA2war1init, model=EWMA2war1mdl, ...){ vars <- exp(pars[1:3]) Q <- model$Q Q[1,1,] <- vars[1] Q[2,2,] <- vars[2] # P1w <- as.matrix(nlme::pdLogChol(pars[4:6])) model$P1[3:4, 3:4] <- P1w Q[3:4, 3:4, beforeWar] <- P1w model$Q <- Q model$H[] <- vars[3] model } EWMA2war1fit <- fitSSM(EWMA2war1mdl, inits=EWMA2war1init, updatefn=EWMA2war1updt) logLiks[7, ] <- logLik2(EWMA2war1fit) logLiks[1:7, ]

This looks like more convergence problems.

Let’s check to confirm we can get logLik(EWMA2fit$model) as before:

{r} EWMA2war1init0 <- c(EWMA2fit$optim.out$par, lnP1.33=-999, lnP1.44=-999, lnP1.34=0) EWMA2war1fit0 <- fitSSM(EWMA2war1mdl, inits=EWMA2war1init0, updatefn=EWMA2war1updt) logLik2(EWMA2war1fit0)

(say something about marginal vs. diffuse likelihood)

Discussion

The following table summarizes the log(likelihood) and the root mean square error or deviation (rmse) from the one-step ahead prediction for each of the nine models described above with a brief comment on each fit:

Model state vector components time varying components parameters to estimate logLik rmse comments
1 EWMA1: Standard Bayesian EWMA 1:level 2:lnQ, lnH 372.9323 0.046 obviously need slope
2 EWMA2: double EWMA 2:level, slope 3:lnQ.lvl, lnQ.slope, lnH 384.6971 0.0427 slope too stiff
3 EWMA2pres1: EWMA2 allowing slope to change only between presidents 2:level, slope Q 3:lnQ.lvl, lnQ.slope, lnH 384.6971 0.0427 same as EWMA2
4 EWMA2pres2: EWMA2 assuming no continuity in slope between presidents 2:level, slope Q, T 3:lnQ.lvl, lnQ.slope, lnH 389.576 0.0421 slightly better than previous models
5 EWMA2pres3: EWMA2 with an adjustment for each president 3:level, slope, Pres Q, T 4:lnQ.lvl, lnQ.slope, lnQ.Pres, lnH 391.4337 0.0411
6 EWMA2war0: EWMA2 with fixed effect of war 4:level, slope, war0, war1 T 3:lnQ.lvl, lnQ.slope, lnH 384.6971 0.0428 convergence problems
7 EWMA2war1: EWMA2 with random effect of war 4:level, slope, war0, war1 Q, T 6:lnQ.lvl, lnQ.slope, lnQ.war00, lnQ.war11, lnQ.war01, lnH ??? ??? more convergence problems

References

[1] D. Bates, M. Mächler, B. M. Bolker, et al. Fitting Linear Mixed-Effects Models using lme4. Tech. rep. Comprehensive R Archive Network (CRAN), 2014.

[2] J. Durbin and S. J. Koopman. Time Series Analysis by State Space Methods, 2nd ed. New York: Oxford U. Pr., 2012.

[3] J. Durbin and S. J. Koopman. Time Series Analysis by State Space Methods, 2nd ed. New York: Oxford U. Pr., 2012. Chap. 2.

[4] S. Graves. “Bayes’ Rule of Information and Monitoring in Manufacturing Integrated Circuits”. In: Bayesian Process Monitoring, Control and Optimization. Ed. by B. M. Colosimo and E. del Castillo. Boca Raton, FL: Chapman and Hall, 2007, pp. 187-244.

[5] S. Graves, S. Bisgaard and M. Kulahci. Designing Bayesian EWMA Monitors Using Gage R & R and Reliability Data. Tech. rep. Productive Systems Engineering, 2002.

[6] A. C. Harvey. Forecasting, Structural Time Series Models and the Kalman Filter. New York: Cambridge U. Pr., 1989.

[7] J. Helske. “KFAS: Exponential Family State Space Models in R”. In: Submitted (2016).

[8] S. Horowitz. Herbert Hoover: Father of the New Deal. Tech. rep. Technical Report, no. 122. Cato Institute, 2011.

[9] R. W. McChesney. Blowing the Roof off the Twenty-First Century. originally appeared in Monthly Review 60/5 (Oct. 2008): 1-19. New York: Monthly Review Pr., 2014. Chap. 6. The U.S. Imperial Triangle and Military Spending.

[10] J. C. Pinheiro and D. M. Bates. Mixed-Effects Models in S and S-PLUS. New York: Springer, 2000, pp. 83-89.

[11] R. Pollin and H. Garrett-Peltier. The U.S. Employment Effects of Military and Domestic Spending Priorities. Tech. rep. University of Massachusetts Department of Economics and Polical Economy Research Institute, 2009.