Stochastic_Process.r

“Generate Stochastic Processes and a realization of Time Series”

“Pedro Valls”

“June 20th, 2021”

# Clear workspace

rm(list=ls())

#set local directory to where data set is
setwd("C:/Users/Pedro/Dropbox/EcoIII2021/Livro/R_CODES")

# Load Packages Functions
load_package<-function(x){
  x<-as.character(match.call()[[2]])
  if (!require(x,character.only=TRUE)){
    install.packages(pkgs=x,repos="http://cran.r-project.org")
    require(x,character.only=TRUE)
  }
}
load_package("readxl")

## Loading required package: readxl

load_package("openxlsx")

## Loading required package: openxlsx

load_package("xts")

## Loading required package: xts

## Loading required package: zoo

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

load_package("fBasics")

## Loading required package: fBasics

## Loading required package: timeDate

## Loading required package: timeSeries

## 
## Attaching package: 'timeSeries'

## The following object is masked from 'package:zoo':
## 
##     time<-

load_package("qrmtools")

## Loading required package: qrmtools

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

## 
## Attaching package: 'qrmtools'

## The following object is masked from 'package:timeSeries':
## 
##     returns

load_package("zoo")
load_package("graphics")
load_package("tseries")

## Loading required package: tseries

load_package("stats")
load_package("moments")

## Loading required package: moments

## 
## Attaching package: 'moments'

## The following objects are masked from 'package:timeDate':
## 
##     kurtosis, skewness

load_package("MASS")

## Loading required package: MASS

load_package("MultiRNG")

## Loading required package: MultiRNG

load_package("tidyverse")

## Loading required package: tidyverse

## Error: package or namespace load failed for 'tidyverse' in loadNamespace(j <- i[[1L]], c(lib.loc, .libPaths()), versionCheck = vI[[j]]):
##  namespace 'jsonlite' 1.7.1 is already loaded, but >= 1.7.2 is required

## Installing package into 'C:/Users/Pedro/OneDrive/Documentos/R/win-library/3.6'
## (as 'lib' is unspecified)

## also installing the dependency 'jsonlite'

## package 'jsonlite' successfully unpacked and MD5 sums checked

## Warning: cannot remove prior installation of package 'jsonlite'

## Warning in file.copy(savedcopy, lib, recursive = TRUE):
## problem copying C:\Users\Pedro\OneDrive\Documentos\R\win-
## library\3.6\00LOCK\jsonlite\libs\x64\jsonlite.dll
## to C:\Users\Pedro\OneDrive\Documentos\R\win-
## library\3.6\jsonlite\libs\x64\jsonlite.dll: Permission denied

## Warning: restored 'jsonlite'

## package 'tidyverse' successfully unpacked and MD5 sums checked
## 
## The downloaded binary packages are in
##  C:\Users\Pedro\AppData\Local\Temp\Rtmp06t1wf\downloaded_packages

## Loading required package: tidyverse

## Error: package or namespace load failed for 'tidyverse' in loadNamespace(j <- i[[1L]], c(lib.loc, .libPaths()), versionCheck = vI[[j]]):
##  namespace 'jsonlite' 1.7.1 is already loaded, but >= 1.7.2 is required

load_package("cowplot")

## Loading required package: cowplot

load_package("urca")

## Loading required package: urca

load_package("vars")

## Loading required package: vars

## Loading required package: strucchange

## Loading required package: sandwich

## Loading required package: lmtest

load_package("tsDyn")

## Loading required package: tsDyn

load_package("ggplot2")

## Loading required package: ggplot2

# Load Library
library(data.table)

## 
## Attaching package: 'data.table'

## The following objects are masked from 'package:xts':
## 
##     first, last

library(ggplot2)
library(tseries)
library(lmtest)
library(Metrics)

## 
## Attaching package: 'Metrics'

## The following object is masked from 'package:tsDyn':
## 
##     mse

library(stargazer)

## 
## Please cite as:

##  Hlavac, Marek (2018). stargazer: Well-Formatted Regression and Summary Statistics Tables.

##  R package version 5.2.2. https://CRAN.R-project.org/package=stargazer

library(xtable)

## 
## Attaching package: 'xtable'

## The following object is masked from 'package:timeSeries':
## 
##     align

## The following object is masked from 'package:timeDate':
## 
##     align

library(forecast)

## 
## Attaching package: 'forecast'

## The following object is masked from 'package:Metrics':
## 
##     accuracy

library(xtable)
library(FinTS)

## 
## Attaching package: 'FinTS'

## The following object is masked from 'package:forecast':
## 
##     Acf

Fix sample sixe

n <-100

Initialize a matrix \(y\) with dimension 100

y <- matrix(data=0, nrow=100, ncol=100)

Initialize a matrix \(z\) with dimension 100

z <- matrix(data=0, nrow=100, ncol=100)

Initialize a matrix \(u\) with dimension 100

u <- matrix(data=0, nrow=100, ncol=100)

Initialize a matrix \(x\) with dimension 100

x <- matrix(data=0, nrow=100, ncol=100)

Initialize the vectors \(x1\), \(x2\), \(x3\) e \(x4\) with the realization of the time series \(x\).

x1 <- matrix(data=0, nrow=100, ncol=1)
x2 <- matrix(data=0, nrow=100, ncol=1)
x3 <- matrix(data=0, nrow=100, ncol=1)
x4 <- matrix(data=0, nrow=100, ncol=1)

Initialize the \(tsx\) vector that will contain the observable time series

tsx <- matrix(data=0, nrow=100, ncol=1)

Fix the random seed for \(y(t,w)\)

set.seed(123456) 

Generate a multivariade normal distribuition with mean = (0,…,0), \(\Sigma\) = diag(1.0,…,1.0)

using mvrnorm para \(y(t,w)\) with order 100

Sigmay <- diag(rep(1,100))
y <- mvrnorm(n, rep(0,100), Sigmay, tol = 1e-6, empirical = TRUE)

Fix random seed for \(z(t,w)\)

set.seed(654321)

Generate a multivariade normal distribuition with mean = (10,…,10), \(\Sigma\) = diag(5.0,…,5.0)

using mvrnorm para \(z(t,w)\) with order 100

Sigmaz <- diag(rep(5,100))
z <- mvrnorm(n, rep(10,100), Sigmaz, tol = 1e-6, empirical = TRUE)

Fix random seed for \(u(t,w)\)

set.seed(123)

Generate the process \(u(t,w)\) using uniform between [\(-\pi\), \(\pi\)]

for(i in 1:n){
  for(j in 1:n){
    u1 <- runif(1, min=-pi,max=pi)
    u[i,j] =u1
  }
} 

Generate \(x(t,w)=y(t,w)*cos(z(t,w)+t*u(t,w))\)

for(i  in  1:n){
  for(j in 1:n){
    x[i,j]=y[i,j]*cos(z[i,j]+i*u[i,j])
  }} 

Select the first four columns are the first four realizations of the stochastic process \(x(t,w)\)

x1 <- ts(x[,1]) 
x2 <- ts(x[,2])
x3 <- ts(x[,3])
x4 <- ts(x[,4])

Transform to data table

stoch.plot <- data.table('id'=1:100, 'x1'= x1, 'x2'=x2, 'x3'=x3, 'x4'=x4)

Plot the four realizations of the time series

ggplot(stoch.plot)+
  geom_line(aes(x=id, y=x1,colour="x1"),linetype="solid") +
  geom_line(aes(x=id, y=x2,colour="x2"),linetype="solid") +
  geom_line(aes(x=id, y=x3,colour="x3"),linetype="solid") +
   geom_line(aes(x=id, y=x4,colour="x4"),linetype="solid") +
  scale_color_manual(values=c("x1"="#FF0000",
                                "x2"="blue",
                              "x3"="black",
                              "x4"="green"
                              ))+
  labs(
  title="Four realizations of the time series x(t,w)",
  x=NULL,
  y=NULL,
  colour = "Legenda"
  ) +
  theme_bw()

## Don't know how to automatically pick scale for object of type ts. Defaulting to continuous.

Generate the \(tsx\) as \(1^{rt}\) obs \(x1\), \(2^{nd}\) obs \(x2\), \(3^{rd}\) obs \(x3\), \(4^{th}\) obs \(x4\),…

teste <- 0
teste1 <- 0
tsx <- 0



for(i  in  1:n){
  teste[i] = floor((i)/4)
  teste1[i]=i-(teste[i])*4
  
  if(teste1[i] == 1){
    tsx[i] <- x1[i] } 
  else if(teste1[i] == 2){
    tsx[i] <- x2[i] } 
  else if(teste1[i] == 3){
    tsx[i] <- x3[i]} 
  else if(teste1[i] == 0){
    tsx[i] <- x4[i]} 
  else {
   
  }
}


tsx = ts(tsx)

Plot of the observable the time series

ggplot(stoch.plot)+
  geom_line(aes(x=id, y=tsx,colour="tsx"),linetype="solid") +
  scale_color_manual(values=c("tsx"="red"
                              ))+
  labs(
    title="Observable time series x(t,w)",
    x=NULL,
    y=NULL,
    colour = "Legenda"
  ) +
  theme_bw()

## Don't know how to automatically pick scale for object of type ts. Defaulting to continuous.