The two stock companies used for the VAR model were tesla and Panasonic. These companies were chosen because they have long term investments with each other for their battery supply and demand needs.

# Lib and data ------------------------------------------------------------

library('tidyverse')
## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --
## v ggplot2 3.3.5     v purrr   0.3.4
## v tibble  3.1.6     v dplyr   1.0.8
## v tidyr   1.2.0     v stringr 1.4.0
## v readr   2.1.2     v forcats 0.5.1
## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library('lubridate')
## 
## Attaching package: 'lubridate'
## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union
library('forecast')
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library('feasts')
## Loading required package: fabletools
## 
## Attaching package: 'fabletools'
## The following objects are masked from 'package:forecast':
## 
##     accuracy, forecast
library('fpp3') 
## -- Attaching packages -------------------------------------------- fpp3 0.4.0 --
## v tsibble     1.1.1     v fable       0.3.1
## v tsibbledata 0.4.0
## -- Conflicts ------------------------------------------------- fpp3_conflicts --
## x lubridate::date()      masks base::date()
## x dplyr::filter()        masks stats::filter()
## x fabletools::forecast() masks forecast::forecast()
## x tsibble::intersect()   masks base::intersect()
## x tsibble::interval()    masks lubridate::interval()
## x dplyr::lag()           masks stats::lag()
## x tsibble::setdiff()     masks base::setdiff()
## x tsibble::union()       masks base::union()
library("seasonal")
## 
## Attaching package: 'seasonal'
## The following object is masked from 'package:tibble':
## 
##     view
library('gridExtra') 
## 
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
## 
##     combine
tsla <- read.csv("TSLA.csv")
pan <- read.csv("PCRFY.csv")

head(tsla)
##         Date  Open  High   Low Close Adj.Close    Volume
## 1 2010-07-01 5.000 5.184 2.996 3.988     3.988 322879000
## 2 2010-08-01 4.100 4.436 3.478 3.896     3.896  75191000
## 3 2010-09-01 3.924 4.632 3.900 4.082     4.082  90229500
## 4 2010-10-01 4.138 4.374 4.000 4.368     4.368  32739000
## 5 2010-11-01 4.388 7.200 4.210 7.066     7.066 141575500
## 6 2010-12-01 7.174 7.284 5.000 5.326     5.326 184464500
head(pan)
##         Date  Open  High   Low Close Adj.Close  Volume
## 1 2010-05-01 14.52 14.70 12.35 12.76  12.69507 6841300
## 2 2010-06-01 12.76 14.06 12.43 12.53  12.46624 6435100
## 3 2010-07-01 12.50 13.55 12.35 13.26  13.19253 6926000
## 4 2010-08-01 13.31 13.35 12.14 12.61  12.54584 6982700
## 5 2010-09-01 12.85 13.80 12.75 13.58  13.51090 8263100
## 6 2010-10-01 13.55 15.00 13.19 14.50  14.42622 5994400
# data exploration and transformation -------------------------------------
tsla$Date <- yearmonth(tsla$Date)
tsla<- as_tsibble(tsla, index = Date)

tsla <- tsla %>% 
  filter(year(Date)>= 2017)

tsla %>% 
  ggplot(aes(x = Date, y = Adj.Close))+
           geom_line()

pan$Date <- yearmonth(pan$Date)

pan <- pan %>% 
  filter(year(Date) >= 2017)

pan <- as_tsibble(pan, index = Date)

pan %>% 
  ggplot(aes(x = Date, y = Adj.Close)) +
  geom_line()

# differnecing ------------------------------------------------------------


pan <- pan %>% 
  mutate(diff.pan = difference(Adj.Close))
pan %>% 
  ggplot(aes(x = Date, y = diff.pan)) +
  geom_line()
## Warning: Removed 1 row(s) containing missing values (geom_path).

tsla <- tsla %>% 
  mutate(diff.tsla = difference(Adj.Close))
tsla %>% 
  ggplot(aes(x = Date, y = diff.tsla))+
  geom_line()
## Warning: Removed 1 row(s) containing missing values (geom_path).

stocks <- bind_cols(tsla, pan$diff.pan) 
## New names:
## * `` -> `...9`
head(stocks) 
## # A tsibble: 6 x 9 [1M]
##       Date  Open  High   Low Close Adj.Close    Volume diff.tsla   ...9
##      <mth> <dbl> <dbl> <dbl> <dbl>     <dbl>     <dbl>     <dbl>  <dbl>
## 1 2017 Jan  43.0  51.7  42.2  50.4      50.4 503398000    NA     NA    
## 2 2017 Feb  50.6  57.5  48.4  50.0      50.0 597700000    -0.388  0.51 
## 3 2017 Mar  50.8  56.4  48.6  55.7      55.7 535176500     5.66   0.42 
## 4 2017 Apr  57.4  63.0  56.9  62.8      62.8 584753000     7.15   0.610
## 5 2017 May  63.0  68.6  58.2  68.2      68.2 740231500     5.39   0.92 
## 6 2017 Jun  68.8  77.4  66.8  72.3      72.3 929755500     4.12   0.710
names(stocks)[names(stocks)== "...9"] <- "diff.pan"
# Var Model  --------------------------------------------------------------

fit <- stocks %>% 
  model(aicc = VAR(vars(diff.tsla, diff.pan)),
        bic = VAR(vars(diff.tsla,diff.pan), ic = "bic"))
report(fit) 
## Warning in report.mdl_df(fit): Model reporting is only supported for individual
## models, so a glance will be shown. To see the report for a specific model, use
## `select()` and `filter()` to identify a single model.
## # A tibble: 2 x 6
##   .model sigma2        log_lik   AIC  AICc   BIC
##   <chr>  <list>          <dbl> <dbl> <dbl> <dbl>
## 1 aicc   <dbl [2 x 2]>   -405.  849.  871.  891.
## 2 bic    <dbl [2 x 2]>   -432.  872.  872.  880.
fit %>%
  augment() %>%
  ACF(.innov) %>%
  autoplot()

fit %>%
  select(aicc) %>%
  forecast() %>%
  autoplot(stocks)
## Warning: Removed 2 row(s) containing missing values (geom_path).

The Panasonic data has some autocrelation which is shown more in the AICc ACF plot. It is still visible in the BIC ACF plot