The two stocks I will use are Chegg (CHGG) and Microsoft Corporation (MSFT). Chegg is an education technology company, and Microsoft is a technology company that sells and supports electronics, software, etc. I assumed that both have seen growth over the past several months as much of education is being delivered via virtual platforms.

library(quantmod)
## Warning: package 'quantmod' was built under R version 4.0.3
## Loading required package: xts
## Warning: package 'xts' was built under R version 4.0.2
## Loading required package: zoo
## Warning: package 'zoo' was built under R version 4.0.2
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## Loading required package: TTR
## Warning: package 'TTR' was built under R version 4.0.2
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
## Version 0.4-0 included new data defaults. See ?getSymbols.
library(fpp2)
## Warning: package 'fpp2' was built under R version 4.0.3
## -- Attaching packages ------------------------------------ fpp2 2.4 --
## v ggplot2   3.3.0     v fma       2.4  
## v forecast  8.12      v expsmooth 2.3
## Warning: package 'forecast' was built under R version 4.0.2
## Warning: package 'fma' was built under R version 4.0.3
## Warning: package 'expsmooth' was built under R version 4.0.3
## 
library(TTR)

Load the data and make it a time series.

CHGG <- read.csv("C:/Users/burtkb/Downloads/CHGG.csv")
MSFT <- read.csv("C:/Users/burtkb/Downloads/MSFT.csv")
chgg.ts = ts(CHGG[,2],start=c(2015,11,30),end=c(2020,11,30),frequency=12)
msft.ts = ts(MSFT[,2],start=c(2015,11,30),end=c(2020,11,30),frequency=12)

#plot the data separately
par(mfrow=c(1,2))
plot.ts(chgg.ts)
plot.ts(msft.ts)

It looks like there is an upward trend and positive correlation which was expected.

#check autocorrelations
ggAcf(chgg.ts)

ggPacf(chgg.ts)

ggAcf(msft.ts)

ggPacf(msft.ts)

#take difference
c.diff = diff(chgg.ts)
m.diff = diff(msft.ts)
ccf(c.diff,m.diff, type = c("correlation","covariance"))

d1 = data.frame(c.diff, m.diff)
d.ts = ts(d1, start=c(2015,11,30),end=c(2020,11,30),frequency=12)
autoplot(d.ts)

d.diff = diff(d.ts)
autoplot(d.diff)

Build VAR Model

library(vars)
## Warning: package 'vars' was built under R version 4.0.3
## Loading required package: MASS
## 
## Attaching package: 'MASS'
## The following objects are masked from 'package:fma':
## 
##     cement, housing, petrol
## Loading required package: strucchange
## Warning: package 'strucchange' was built under R version 4.0.3
## Loading required package: sandwich
## Warning: package 'sandwich' was built under R version 4.0.2
## Loading required package: urca
## Warning: package 'urca' was built under R version 4.0.2
## Loading required package: lmtest
## Warning: package 'lmtest' was built under R version 4.0.2
VARselect(d.diff,lag.max=8,type="const")[["selection"]]
## AIC(n)  HQ(n)  SC(n) FPE(n) 
##      8      8      8      8
var1 = VAR(d.diff, p=1, type="const")
summary(var1)
## 
## VAR Estimation Results:
## ========================= 
## Endogenous variables: c.diff, m.diff 
## Deterministic variables: const 
## Sample size: 59 
## Log Likelihood: -415.26 
## Roots of the characteristic polynomial:
## 0.9302 0.6526
## Call:
## VAR(y = d.diff, p = 1, type = "const")
## 
## 
## Estimation results for equation c.diff: 
## ======================================= 
## c.diff = c.diff.l1 + m.diff.l1 + const 
## 
##           Estimate Std. Error t value Pr(>|t|)  
## c.diff.l1  -0.1719     0.3943  -0.436   0.6645  
## m.diff.l1  -0.3431     0.1683  -2.039   0.0462 *
## const      -1.2643     1.4055  -0.900   0.3722  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 10.71 on 56 degrees of freedom
## Multiple R-Squared: 0.4918,  Adjusted R-squared: 0.4736 
## F-statistic: 27.09 on 2 and 56 DF,  p-value: 5.89e-09 
## 
## 
## Estimation results for equation m.diff: 
## ======================================= 
## m.diff = c.diff.l1 + m.diff.l1 + const 
## 
##           Estimate Std. Error t value Pr(>|t|)    
## c.diff.l1   1.0623     0.8934   1.189 0.239418    
## m.diff.l1  -1.4109     0.3813  -3.700 0.000493 ***
## const      -2.9061     3.1849  -0.912 0.365436    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 24.28 on 56 degrees of freedom
## Multiple R-Squared: 0.5224,  Adjusted R-squared: 0.5053 
## F-statistic: 30.62 on 2 and 56 DF,  p-value: 1.035e-09 
## 
## 
## 
## Covariance matrix of residuals:
##        c.diff m.diff
## c.diff  114.8  250.4
## m.diff  250.4  589.4
## 
## Correlation matrix of residuals:
##        c.diff m.diff
## c.diff 1.0000 0.9628
## m.diff 0.9628 1.0000
serial.test(var1, type="PT.asymptotic")
## 
##  Portmanteau Test (asymptotic)
## 
## data:  Residuals of VAR object var1
## Chi-squared = 64.488, df = 60, p-value = 0.3226
fv1 = forecast(var1)
autoplot(fv1)

tsdisplay(residuals(var1))