Introduction

Data Exploration

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   72.21   94.86  103.67  104.26  114.20  134.75

Work from Midterm Project

## 
## Call:
## arima(x = IPI, order = c(0, 1, 1), seasonal = list(order = c(0, 1, 1), period = 12))
## 
## Coefficients:
##           ma1     sma1
##       -0.2907  -1.0000
## s.e.   0.0690   0.1252
## 
## sigma^2 estimated as 21.39:  log likelihood = -616.27,  aic = 1238.53

POMP Model Description

Background

  • There are some observations of Gross Domestic Product (GDP) being a GBM. (Conover, n.d.) Since IPI is closely related to GDP, it is reasonable to use GBM to model IPI.

  • A GBM \(\{S_t\}\) is an exponentiated Brownian motion defined through \[\frac{dS_t}{S_t} = \mu dt + \sigma dW_t\] where \(W_t\) is a Brownian motion. (Avdiu, n.d.)

  • There is a recursive procedure for simulating values of \(\{S_t\}\): \[S_{t+1} = S_t \exp( (\mu-\frac{1}{2}\sigma^2) (t_{i+1} - t_i) + \sigma \sqrt{t_{i+1} - t_i} Z_i)\] with \(0 < t_0 < t_1 < ...\) and \(Z_i \sim N(0,1)\). (Avdiu, n.d.)

Our Model

  • In our monthly data, \(t_{i+1} - t_i = 1\).

  • We present a POMP model with GBM as \[S_{n+1} = S_n \exp((\mu-\frac{1}{2}\sigma^2) + \sigma E_n)\] where \(S\) and \(E\) are state variables, and \(E_n \sim N(0, \phi^2)\).

  • The measurement model is \[IPI_n | S_n \sim N(0,S_n)\].

  • We initialize \(S_0 = 1\) and \(E_0 = 0\).

  • The following is the corresponding pomp implementation.

ipi_statenames <- c("S","E")
ipi_paramnames <- c("mu","sigma","phi")
ipi_obsnames <- "IPI"

ipi_step <- Csnippet("
E = rnorm(0,phi);
S = S*exp((mu-sigma*sigma/2)+sigma*E);
")

ipi_init <- Csnippet("
S=1;
E=0;
")

ipi_fromEstimationScale <- "
 Tmu = exp(mu);
 Tsigma = exp(sigma);
 Tphi = exp(phi);
"

ipi_toEstimationScale <- "
 Tmu = log(mu);
 Tsigma = log(sigma);
 Tphi = log(phi);
"

dmeas <- Csnippet("lik = dnorm(IPI,0,S,give_log);")
rmeas <- Csnippet("IPI = rnorm(0,S);")

pomp(dt,times="Date",t0=0,
     rprocess=discrete.time.sim(step.fun=ipi_step,delta.t=1),
     initializer=ipi_init,rmeasure=rmeas,dmeasure=dmeas,
     statenames=ipi_statenames,
     paramnames=ipi_paramnames,
     fromEstimationScale=Csnippet(ipi_fromEstimationScale),
    toEstimationScale=Csnippet(ipi_toEstimationScale)
     ) -> ipi_pomp

Run Level

  • We use three run levels to develop and refine the model.
run_level <- 3
switch(run_level,
       {ipi_Np=100; ipi_Nmif=10; ipi_Neval=10; ipi_Nglobal=10; ipi_Nlocal=10}, 
       {ipi_Np=20000; ipi_Nmif=100; ipi_Neval=10; ipi_Nglobal=10; ipi_Nlocal=10}, 
       {ipi_Np=60000; ipi_Nmif=300; ipi_Neval=10; ipi_Nglobal=100; ipi_Nlocal=20}
)

Likelihood Slice

  • To have an idea of what the likelihood surface looks like, we design the following slice.
sliceDesign(
  c(mu=0.1,sigma=0.2,phi=0.5),
  mu=rep(seq(from=-0.5,to=0.5,length=50),each=3),
  sigma=rep(seq(from=0,to=4,length=50),each=3),
  phi=rep(seq(from=0,to=4,length=50),each=3)
  ) -> p

  • Along \(\mu\) direction, the log likelihood is maximized when \(\mu\) is about 0.1.

  • Along \(\sigma\) direction, the log likelihood is maximized when \(\sigma\) is about 0.5.

  • Along \(\phi\) direction, the log likelihood is maximized when \(\phi\) is about 2.

Running a Partical Filter

  • Using \(\mu=0.1\), \(\sigma=0.5\) and \(\phi=2\) as point estimate, we run a basic particle filter.
##                          se 
## -1.453785e+03  8.526962e-02
  • The log likelihood estimate is about -1453.78 with a Monte Carlo standard error of 0.08.

Diagnostic

  • The fact that both local search and global search give similar max log likelihood shows confidence in our maximization procedure.

  • We also look at the diagnostic plots.

plot(mifs_global)

  • The diagnostic plots show that the parameter estimations are quite stable.

Conclusion

Acknowledgment

Reference

[1] Avdiu, K. Brownian Motion & Geometric Brownian Motion. Retrieved April 26, 2018, from http://homepage.univie.ac.at/kujtim.avdiu/dateien/BrownianMotion.pdf.

[2] Conover, J. Historical Economics. Retrieved April 26, 2018, from http://www.johncon.com/john/historical.economics/.

[3] Huerta, G., Lopes, H. F. 1999. Bayesian forecasting and inference in latent structure for the Brazilian GDP and Industrial Production Index. Retrieved April 26, 2018, from ftp://ftp.stat.duke.edu/pub/WorkingPapers/99-08.html.

[4] Industrial Production Index. Retrieved April 26, 2018, from https://fred.stlouisfed.org/series/INDPRO.

[5] Ionides, E. Stats 531 (Winter 2018) ‘Analysis of Time Series’. Retrieved April 26, 2018, from https://ionides.github.io/531w18/.

[6] West, M. SOME TIME SERIES DATA SETS. Retrieved April 26, 2018, from http://www2.stat.duke.edu/~mw/ts_data_sets.html.