---
title: "VisitorCounts"
output: 
  rmarkdown::html_vignette:
    toc: true
vignette: >
  %\VignetteIndexEntry{national_park_analysis}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

In this vignette, functions in the VisitorCounts package are demonstrated using park visitation data from Yellowstone National Park.

# Sample Datasets: `park_visitation` and `flickr_userdays`

First, we load two datasets: `park_visitation` stores 156 monthly observations spanning 2005 through 2017 of flickr user-days (PUD) and visitor counts by the national park service (NPS) for 20 popular national parks in the United States. Second, `flickr_userdays` stores log US flickr user-days for the corresponding time period.


```{r}
library(VisitorCounts)
data("park_visitation")
data("flickr_userdays")
```

## Sample data for Yellowstone National Park

For the purposes of this vignette, three time series are extracted from these datasets. First, `log_yellowstone_pud` is a time series of 156 monthly observations of flickr photo-user-days geolocated within Yellowstone National Park. Second, `log_yellowstone_nps` is a time series of 156 monthly observations of counts of park visitation by the national park service. Third, `flickr_userdays` is a time series of 156 monthly observations of log flickr user-days taken within the United States.  



```{r, fig.width = 7, fig.height = 5}



yellowstone_pud <- park_visitation[park_visitation$park == "YELL",]$pud #photo user days
yellowstone_nps <- park_visitation[park_visitation$park == "YELL",]$nps #national park service counts

yellowstone_pud <- ts(yellowstone_pud, start = 2005, freq = 12)
yellowstone_nps <- ts(yellowstone_nps, start = 2005, freq = 12)

log_yellowstone_pud <- log(yellowstone_pud)
log_yellowstone_nps <- log(yellowstone_nps)

log_flickr_userdays <- log(flickr_userdays)
```


```{r, fig.width = 7, fig.height = 5}
plot(log_yellowstone_pud, main = "Yellowstone Photo-User-Days (PUD)", ylab = "PUD")
plot(log_yellowstone_nps, main = "Yellowstone National Park Service Visitation Counts (NPS)", ylab = "NPS")
plot(log_flickr_userdays, main = "Log US Flickr user-days", ylab = "UD")
```




# visitation_model()

The `visitation_model()` function uses social media data, such as the log flickr photo-user-days in `log_yellowstone_pud`, coupled with a popularity measure of the social media platform, like the log US flickr userdays in `log_flickr_userdays`, to model percent changes in visitation counts.
By default, `visitation_model()` assumes that no visitation counts are available. 
```{r}
yell_visitation_model <- visitation_model(log_yellowstone_pud,
                                          log_flickr_userdays, is_output_logged = TRUE, is_input_logged = TRUE)
```
If national park data is available, a reference series may be supplied to assist in parameter estimates:

```{r}
yell_visitation_model_nps <- visitation_model(log_yellowstone_pud,
                                              log_flickr_userdays,
                                              ref_series = log_yellowstone_nps, is_output_logged = TRUE, is_input_logged = TRUE)
```



## plot.visitation_model()

By default, `plot.visiation_model()` plots the differenced series. Typical graphical parameters may be passed to `plot.visitation_model()`, such as line width:

```{r, fig.width = 7, fig.height = 5}
true_differences <- diff(log_yellowstone_nps)
lower_bound <- min(c(true_differences,diff(yell_visitation_model$visitation_fit)))-1
upper_bound <- max(c(true_differences,diff(yell_visitation_model$visitation_fit)))

plot(yell_visitation_model, ylim = c(lower_bound, upper_bound), lwd = 2)
lines(diff(log_yellowstone_nps), col = "red")
legend("bottom",c("Model Fit","True Differences"),col = c("black","red"),lty = c(1,1))
```

```{r, fig.width = 7, fig.height = 5}
true_differences <- diff(log_yellowstone_nps)
lower_bound <- min(c(true_differences,diff(yell_visitation_model_nps$visitation_fit)))-1
upper_bound <- max(c(true_differences,diff(yell_visitation_model_nps$visitation_fit)))

plot(yell_visitation_model_nps, ylim = c(lower_bound, upper_bound), 
     lwd = 2,
     main = "Fitted Values for Visitation Model (NPS assisted)", difference = TRUE)
lines(true_differences, col = "red")
legend("bottom",c("Model Fit","True Differences"),col = c("black","red"),lty = c(1,1))
```

## summary.visitation_model()

Parameters can be inspected using `summary.visitation_model()`. Two examples can be seen below:

```{r}
summary(yell_visitation_model)
summary(yell_visitation_model_nps)
```


## predict.visitation_model()

Forecasts can be made using `predict.visitation_model()`, whose output is a `visitation_forecast` class object which can be inspected using `plot` or `summary` functions.

```{r}
yellowstone_visitation_forecasts <- predict(yell_visitation_model, n_ahead = 12)
yellowstone_visitation_forecasts_nps <- predict(yell_visitation_model_nps, n_ahead = 12)

yellowstone_visitation_forecasts_withpast <- predict(yell_visitation_model, n_ahead = 12, only_new = FALSE)
```


### plot.visitation_forecast()

Forecasts can be plotted using `plot.visitation_forecast()`:

```{r, fig.width = 7, fig.height = 5}
plot(yellowstone_visitation_forecasts, difference = TRUE)
plot(yellowstone_visitation_forecasts_nps, main = "Forecasts for Visitation Model (NPS Assisted)", date_label = "%b", date_breaks = "1 month")


plot(yellowstone_visitation_forecasts_withpast, difference = TRUE, date_breaks = "1 year", date_label = "%y")
```



### summary.visitation_forecast()

```{r}
summary(yellowstone_visitation_forecasts)
summary(yellowstone_visitation_forecasts_nps)
```


# auto_decompose()

The automatic decomposition function uses singular-spectrum analysis, as implemented by the Rssa package, in conjunction with an automated procedure for classifying components to decompose a time series into trend, seasonality and noise. 

```{r}
yell_pud_decomposition <- auto_decompose(yellowstone_pud)
```


## plot.decomposition()

Several plot options are available for examining this decomposition.

```{r, fig.width = 7, fig.height = 5}
plot(yell_pud_decomposition)

plot(yell_pud_decomposition, type = "period")

plot(yell_pud_decomposition, type = "classical")
```


## summary.decomposition()

The eigenvector grouping can be examined using `summary.decomposition`.

```{r}
summary(yell_pud_decomposition)
```

## predict.decomposition()

Forecasts can be made using `predict.decomposition()`:

```{r, fig.width = 7, fig.height =5}
plot(predict(yell_pud_decomposition, n_ahead = 12)$forecast, main = "Decomposition 12-ahead Forecast", ylab = "Forecast Value")
```