---
title: "Time series with a sampling error component"
author: "Jean Palate"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Vignette Title}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

## Introduction

This document reproduces with rjd3sts the example N°6 of the paper:

*Bell W.R. (2011). REGCMPNT - A Fortran Program for Regression Models with ARIMA Component Errors.
Journal of Statistical Software. Volume 41, Issue 7.*

Details on the considered model can be retrieved from the original paper (pages 17-20).



```{r echo=FALSE}
suppressPackageStartupMessages(library(rjd3toolkit))
suppressPackageStartupMessages(library(rjd3sts))

# Data

nr055<-c(115.2,110.4,111.5,127.9,150,149.5,139.5,144.6,176,173.2,172.7,143.5,135.3,174.5,164.5,179.1,
         164.8,168.7,115.8,163.7,156.3,129.3,125.5,120.8,84.4,93.6,83,120.5,110.9,115.3,108.5,107.7,
         99.8,98.7,107.7,102.1,95.9,92.2,99.3,92.8,127.4,105.2,108,118.1,108.5,133.2,127.3,110.1,103.1,
         99.8,101.1,105.8,101.6,110.9,122.9,129.8,129.4,149.7,141.4,125,122.3,110.3,108.9,122.8,115.8,
         116.9,134.2,158.8,141.9,174.1,178.6,154.5)

h<-c(.042,.042,.067,.122,.129,.135,
     .152,.168,.173,.177,.179,.179,
     .179,.182,.179,.177,.175,.170,
     .162,.165,.152,.149,.149,.135,
     .144,.140,.144,.159,.152,.149,
     .144,.129,.144,.140,.140,.149,
     .149,.156,.152,.152,.156,.165,
     .159,.159,.159,.149,.144,.144,
     .152,.159,.156,.156,.162,.156,
     .165,.162,.165,.168,.170,.159,
     .144,.149,.144,.144,.129,.135,
     .140,.135,.149,.140,.140,.115)

```


### Definition of the model in rjd3sts

Remark: the code presented below will be further simplified (the definition of the equation will be suppressed, as it is already the case for non-weighted state components).

```{r}

y<-log(nr055)

# create the model
model<-rjd3sts::model()
# create the components and add them to the model

# airline block
airline<-rjd3sts::sarima("airline", 12, c(0,1,1), c(0,1,1))
# fixed AR model for the sampling error
ar<-rjd3sts::ar("ar", ar=c(0.6, 0.246), fixedar=TRUE, variance=0.34488, fixedvariance=TRUE)

rjd3sts::add(model, airline)
rjd3sts::add(model, ar)

# creation of the unique equation
eq<-rjd3sts::equation("eq")
rjd3sts::add_equation(eq, "airline")
rjd3sts::add_equation(eq, "ar", loading = rjd3sts::var_loading(pos=0, weights = h))
rjd3sts::add(model, eq)

# estimate the model
# it is important to note that we don't use the concentrated likelihood (which is the default option,
# faster and more stable)
# In this case the optimization procedure is BFGS instead of Levenberg-Marquardt.

rslt<-rjd3sts::estimate(model, y, concentrated=FALSE)


```

### Main results

```{r main, echo=FALSE}
p<-rjd3toolkit::result(rslt, 'parameters')
q<-rjd3sts::smoothed_components(rslt)
qe<-rjd3sts::smoothed_components_stdev(rslt)

```

#### Estimated parameters

```{r parameter, echo=FALSE}

cat("Airline:\n\n")
cat("Innovation variance: ", p[1], "\n")
cat("theta: ", p[2], "\n")
cat("btheta: ", p[3], "\n")

```
#### Original series (solid) and signal extraction estimates (dashed)
 
```{r signal,fig.width=7, fig.height=4}

data<-cbind(nr055, exp(q[,1]))
matplot(data, type='l')

```

#### Signal extraction estimates of the sampling error
 
```{r samplingerror,fig.width=7, fig.height=4}
noise<-exp(q[,2])
plot(noise, type='l')

```

#### Sampling error and signal extraction error
 
```{r stderr,fig.width=7, fig.height=4}
data<-cbind(h, qe[,2])
matplot(data, type='l')

```

#### CV improvvement due to signal extraction
 
```{r improvement,fig.width=7, fig.height=4}
improvement<-(h-qe[,2])/h*100
plot(improvement, type='l')

```