Simulation and estimation of an integer-valued trawl process

##Trawl processes The trawl package can be used to simulate and estimate univariate and bivariate integer-valued trawl processes as described in (Veraart 2019).

###Simulation and inference in the univariate case A univariate trawl process can be simulated using the function sim_UnivariateTrawl. Currently, one can choose amongst two marginal distributions (Poisson and negative binomial) and four trawl functions (exponential, double-exponential, supIG and long memory).

The simulation of the univariate trawl process follows the algorithm described in Section 4 in (Veraart 2019).

Consider a trawl process Yt = L(At) = X0, t + Xt, where L is a Levy basis, A a trawl. Also, X0, t = L({(x, s) : s ≤ 0, 0 ≤ x ≤ g(s − t)}) and Xt = L({(x, s) : 0 < s ≤ t, 0 ≤ x ≤ g(s − t)}), for a trawl function g. Since X0, t is asymptotically negligible in the sense that it converges to zero as t → ∞, the simulation algorithm focusses on the term Xt only (and introduces a burn-in period).

Note that a realisation of ${\bf L}$ consists of a countable set R of points $({\bf y},x,s)$ in $\mathbb{N}_0^n\setminus\{ {\bf 0}\} \times [0,1]\times \mathbb{R}$. When we project the point pattern to the time axis, we obtain the arrival times of a Poisson process Nt with intensity v = ν(ℝn). The corresponding arrival times are denoted by t1, …, tNt and we associate uniform heights U1, …, UNt with them, see (Barndorff-Nielsen et al. 2014) for a detailed discussion in the univariate case. So as soon as we have specified the jump size distribution of the ${\bf C}$, we can use the representation
$$ X_t= \sum_{j=1}^{N_t}C_j\mathbb{I}_{\{U_j \leq g(t_j-t)\}}, $$ to simulate each component.

We want to simulate X on a Δ-grid of [0, t], where Δ > 0, i.e., we want to find (X0, X1Δ, …, X⌊t/Δ⌋Δ). We proceed as follows:

For the inference, we follow the (generalised) method of moments described in Section 4 (Veraart 2019).

The relevant functions which are available in the trawl package are fit_Exptrawl, fit_DExptrawl fit_supIGtrawl, fit_LMtrawl, for fitting an exponential, double-exponential, supIG or long memory trawl. After the Lebesgue measure of the trawl has been estimated using one of these four functions, the marginal Poisson or negative binomial law can be estimated using fit_marginalPoisson and fit_marginalNB, respectively.

The following examples illustrate the use of these functions.

####Trawl processes with Poisson marginal law

##Poisson with Exp trawl
set.seed(1)
t<-1000
Delta<-1
v<-250
lambda <-0.25

#Simulate a univariate trawl process with exponential trawl function and Poisson marginal law
trawl <-trawl::sim_UnivariateTrawl(t,Delta,burnin=50,marginal =c("Poi"),trawl ="Exp",v=v, lambda1=lambda)

#Plot the sample path of the simulated process
plot(trawl,type="l")


#Fit the exponential trawl function to the simulated data
fittrawlfct <- trawl::fit_Exptrawl(trawl,Delta, plotacf=TRUE,lags=500)


#Fit the marginal Poisson law
fitmarginallaw <-trawl::fit_marginalPoisson(trawl, fittrawlfct$LM, plotdiag=TRUE)


###Print the results
print(paste("lambda: estimated:", fittrawlfct$lambda, ", theoretical:", lambda))
#> [1] "lambda: estimated: 0.241963140792181 , theoretical: 0.25"
print(paste("v: estimated:", fitmarginallaw$v, ", theoretical:", v))
#> [1] "v: estimated: 242.197119149327 , theoretical: 250"
##Poisson with supIG trawl
set.seed(1)
t<-1000
Delta<-1
v<-250
delta <-0.2
gamma <-0.5

#Simulate a univariate trawl process with supIG trawl function and Poisson marginal law
trawl <-trawl::sim_UnivariateTrawl(t,Delta,burnin=100,marginal =c("Poi"),trawl ="supIG",v=v, delta=delta, gamma=gamma)

#Plot the sample path of the simulated process
plot(trawl,type="l")


#Fit the supIG trawl function to the simulated data
fittrawlfct <- trawl::fit_supIGtrawl(trawl,Delta, plotacf=TRUE,lags=500)


#Fit the marginal Poisson law
fitmarginallaw <-trawl::fit_marginalPoisson(trawl, fittrawlfct$LM, plotdiag=TRUE)


###Print the results
print(paste("delta: estimated:", fittrawlfct$delta, ", theoretical:", delta))
#> [1] "delta: estimated: 0.189796214920193 , theoretical: 0.2"
print(paste("gamma: estimated:", fittrawlfct$gamma, ", theoretical:", gamma))
#> [1] "gamma: estimated: 0.511685065248945 , theoretical: 0.5"
print(paste("v: estimated:", fitmarginallaw$v, ", theoretical:", v))
#> [1] "v: estimated: 230.833355301605 , theoretical: 250"

####Trawl processes with negative binomial marginal law

##Negative binomial with Exp trawl
set.seed(1)
t<-1000
Delta<-1
m<-200
theta<-0.5
lambda <-0.25

m*abs(log(1-theta)) #=v
#> [1] 138.6294

#Simulate a univariate trawl process with exponential trawl function and negative binomial marginal law
trawl <-trawl::sim_UnivariateTrawl(t,Delta,burnin=50,marginal =c("NegBin"),trawl ="Exp",m=m, theta=theta, lambda1=lambda)

#Plot the sample path of the simulated process
plot(trawl,type="l")


#Fit the exponential trawl function to the simulated data
fittrawlfct <- trawl::fit_Exptrawl(trawl,Delta, plotacf=TRUE,lags=500)


#Fit the marginal negative binomial law
fitmarginallaw <-trawl::fit_marginalNB(trawl, fittrawlfct$LM, plotdiag=TRUE)


###Print the results
print(paste("lambda: estimated:", fittrawlfct$lambda, ", theoretical:", lambda))
#> [1] "lambda: estimated: 0.281215255331676 , theoretical: 0.25"
print(paste("m: estimated:", fitmarginallaw$m, ", theoretical:", m))
#> [1] "m: estimated: 231.550570999856 , theoretical: 200"
print(paste("theta: estimated:", fitmarginallaw$theta, ", theoretical:", theta))
#> [1] "theta: estimated: 0.490887447419009 , theoretical: 0.5"
##Negative binomial with LM trawl
set.seed(1)
t<-1000
Delta<-1
m<-200
theta<-0.5
m*abs(log(1-theta)) #=v
#> [1] 138.6294
alpha <-2
H <-3

#Simulate a univariate trawl process with long memory trawl function and negative binomial marginal law
trawl <-trawl::sim_UnivariateTrawl(t,Delta,burnin=100,marginal =c("NegBin"),trawl ="LM",m=m,theta=theta, alpha=alpha, H=H)

#Plot the sample path of the simulated process
plot(trawl,type="l")


#Fit the long memory trawl function to the simulated data
fittrawlfct <- trawl::fit_LMtrawl(trawl,Delta, plotacf=TRUE,lags=500)


#Fit the marginal negative binomial law
fitmarginallaw <-trawl::fit_marginalNB(trawl, fittrawlfct$LM, plotdiag=TRUE)


###Print the results
print(paste("alpha: estimated:", fittrawlfct$alpha, ", theoretical:", alpha))
#> [1] "alpha: estimated: 1.87237962845072 , theoretical: 2"
print(paste("H: estimated:", fittrawlfct$H, ", theoretical:", H))
#> [1] "H: estimated: 3.05651225152775 , theoretical: 3"
print(paste("m: estimated:", fitmarginallaw$m, ", theoretical:", m))
#> [1] "m: estimated: 197.982721733593 , theoretical: 200"
print(paste("theta: estimated:", fitmarginallaw$theta, ", theoretical:", theta))
#> [1] "theta: estimated: 0.522796369302733 , theoretical: 0.5"

###Simulation and inference in the bivariate case The simulation and inference methods described above can be generalised to a multivariate setting as described in (Veraart 2019). Currently, the routines are implemented for a bivariate setting. A bivariate trawl process can be simulated using the function sim_BivariateTrawl.

In Section 3 of (Veraart 2019) two types of dependent bivariate trawl processes are discussed: the case of dependence through a common factor (see Example 8 in the underlying paper (Veraart 2019)), which is referred to as fullydep in the options of sim_BivariateTrawl, and the case of dependence through a common factor with additional independent components (see Example 9 in the underlying paper (Veraart 2019)), which is referred to as dep in the options of sim_BivariateTrawl.

The current implementation allows for simulating a bivariate trawl process with either Poisson or negative binomial marginal laws (the mixed case is possible) and the four trawl functions used in the univariate case: exponential, double-exponential, supIG and long memory.

When it comes to the inference, we provide the following functions in addition to the ones for the univariate estimation: The function fit_trawl_intersection can be used to compute the Lebesgue measure of the intersection of two trawls, where the trawl functions are either exponential, double-exponential, supIG or long-memory. More precisely, consider trawls A1 and A2 characterised by trawl functions g1 and g2. Then the function fit_trawl_intersection computes R12(0) := Leb(A1 ∩ A2), where Leb denotes the Lebesgue measure. For the special cases that both trawls are either exponential or long memory, the special functions fit_trawl_intersection_Exp and fit_trawl_intersection_LM are also available.

####Example We illustrate the simulation and estimation of a bivariate trawl process in the following.

#Exponential trawls
lambda1 <- 1.2
lambda2 <- 1.5

#Parameters of the negative binomial marginal law
m1<- 2.1
theta1 <- 0.9
a1 <-27.3

m2 <- 2.3
theta2 <-0.9
a2 <- 35.3


kappa12 <-m1
kappa1 <- 0
kappa2 <- m2-kappa12

#Specify the time period and grid
t <-720
Delta <- 1


##Fix the seed
set.seed(1)

#Simulate the bivariate trawl process with common factor and independent components ("dep") and
#negative binomial marginal law. Both trawl functions are chosen as exponentials.
simdata<-trawl::sim_BivariateTrawl(t, Delta, burnin=10,marginal ="NegBin", dependencetype="dep",
                              trawl1 ="Exp", trawl2 ="Exp",
                              kappa1=kappa1,kappa2=kappa2,kappa12=kappa12,a1=a1,a2=a2,
                              lambda11=lambda1, lambda21 =lambda2)
  
trawl1 <- simdata[,1]
trawl2 <- simdata[,2]

#####Produce a bivariate histogram of the simulated data
trawl::plot_2and1hist(trawl1,trawl2)


#####Produce a bivariate histogram of the simulated data
# using ggplot2
trawl::plot_2and1hist_gg(trawl1,trawl2)
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.



###Fit the parameters of the exponential trawl functions
fit1 <- trawl::fit_Exptrawl(trawl1)
fit2 <- trawl::fit_Exptrawl(trawl2)

###Fit negative binomial model
fitNB1 <- trawl::fit_marginalNB(trawl1,fit1$LM,TRUE)

fitNB2 <- trawl::fit_marginalNB(trawl2,fit2$LM,TRUE)


###Compute the intersection between the two trawls
R12 <- trawl::fit_trawl_intersection("Exp","Exp",lambda11=fit1$memory,lambda21=fit2$memory, LM1=fit1$LM, LM2=fit2$LM,TRUE)
#> [1] "There are no roots"

###Fit the remaining parameters from the bivariate negative binomial law
kappa12_est <- min(stats::cov(trawl1,trawl2)/(R12*fitNB1$a*fitNB2$a),fitNB1$m,fitNB2$m)
kappa1_est <-max(fitNB1$m -kappa12,0)
kappa2_est <-max(fitNB2$m -kappa12,0)

###Print the results
print("Estimated parameters of the exponential trawl functions:")
#> [1] "Estimated parameters of the exponential trawl functions:"
print(paste("lambda1: estimated:", fit1$lambda, ", theoretical:", lambda1))
#> [1] "lambda1: estimated: 1.04367205874298 , theoretical: 1.2"
print(paste("lambda2: estimated:", fit2$lambda, ", theoretical:", lambda2))
#> [1] "lambda2: estimated: 1.37334279532044 , theoretical: 1.5"

print("Estimated parameters of the bivariate negative binomial law:")
#> [1] "Estimated parameters of the bivariate negative binomial law:"
print(paste("m1: estimated:", fitNB1$m, ", theoretical:", m1))
#> [1] "m1: estimated: 1.86651625219138 , theoretical: 2.1"
print(paste("theta1: estimated:", fitNB1$theta, ", theoretical:", theta1))
#> [1] "theta1: estimated: 0.964798547340514 , theoretical: 0.9"
print(paste("alpha1: estimated:", fitNB1$a, ", theoretical:", a1))
#> [1] "alpha1: estimated: 27.4079185502173 , theoretical: 27.3"


print(paste("m2: estimated:", fitNB2$m, ", theoretical:", m2))
#> [1] "m2: estimated: 2.22130914996701 , theoretical: 2.3"
print(paste("theta2: estimated:", fitNB2$theta, ", theoretical:", theta2))
#> [1] "theta2: estimated: 0.972054470109286 , theoretical: 0.9"
print(paste("alpha2: estimated:", fitNB2$a, ", theoretical:", a2))
#> [1] "alpha2: estimated: 34.7838983161416 , theoretical: 35.3"

print(paste("kappa12: estimated:", kappa12_est, ", theoretical:", kappa12))
#> [1] "kappa12: estimated: 1.86651625219138 , theoretical: 2.1"
print(paste("kappa1: estimated:", kappa1_est, ", theoretical:", kappa1))
#> [1] "kappa1: estimated: 0 , theoretical: 0"
print(paste("kappa2: estimated:", kappa2_est, ", theoretical:", kappa2))
#> [1] "kappa2: estimated: 0.121309149967011 , theoretical: 0.2"

##References

Barndorff-Nielsen, Ole Eiler, Asger Lunde, Neil Shephard, and Almut Elisabeth Dorothea Veraart. 2014. “Integer-Valued Trawl Processes: A Class of Stationary Infinitely Divisible Processes.” Scandinavian Journal of Statistics 41 (3): 693–724.
Veraart, Almut Elisabeth Dorothea. 2019. “Modelling, Simulation and Inference for Multivariate Time Series of Counts Using Trawl Processes.” Journal of Multivariate Analysis 169: 110–29. https://doi.org/10.1016/j.jmva.2018.08.012.