## ----setup, include = FALSE--------------------------------------------------- knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) ## ----------------------------------------------------------------------------- library(emhawkes) ## ----------------------------------------------------------------------------- set.seed(1107) mu1 <- 0.3; alpha1 <- 1.2; beta1 <- 1.5 hspec1 <- new("hspec", mu = mu1, alpha = alpha1, beta = beta1) show(hspec1) ## ---- warning=FALSE----------------------------------------------------------- res1 <- hsim(hspec1, size = 1000) summary(res1) ## ----------------------------------------------------------------------------- # first and third columns are the same cbind(res1$lambda[1:5], res1$lambda_component[1:5], mu1 + res1$lambda_component[1:5]) ## ----------------------------------------------------------------------------- # second and third columns are the same cbind(res1$lambda[1:5], res1$rambda[1:5], res1$lambda[1:5] + alpha1) ## ----------------------------------------------------------------------------- # By definition, the following two are equal: res1$lambda[2:6] mu1 + (res1$rambda[1:5] - mu1) * exp(-beta1 * res1$inter_arrival[2:6]) ## ---- warning=FALSE----------------------------------------------------------- logLik(hspec1, inter_arrival = res1$inter_arrival) ## ---- warning=FALSE----------------------------------------------------------- # initial value for numerical optimization mu0 <- 0.5; alpha0 <- 1.0; beta0 <- 1.8 hspec0 <- new("hspec", mu = mu0, alpha = alpha0, beta = beta0) # the intial values are provided through hspec mle <- hfit(hspec0, inter_arrival = res1$inter_arrival) summary(mle) ## ----------------------------------------------------------------------------- mu2 <- matrix(c(0.2), nrow = 2) alpha2 <- matrix(c(0.5, 0.9, 0.9, 0.5), nrow = 2, byrow = TRUE) beta2 <- matrix(c(2.25, 2.25, 2.25, 2.25), nrow = 2, byrow = TRUE) hspec2 <- new("hspec", mu=mu2, alpha=alpha2, beta=beta2) print(hspec2) ## ----------------------------------------------------------------------------- res2 <- hsim(hspec2, size=1000) summary(res2) ## ----------------------------------------------------------------------------- # Under bi-variate model, there are two types, 1 or 2. res2$type[1:10] ## ----------------------------------------------------------------------------- res2$N[1:3, ] ## ----------------------------------------------------------------------------- res2$lambda[1:3, ] ## ----------------------------------------------------------------------------- res2$lambda_component[1:3, ] ## ----------------------------------------------------------------------------- mu2[1] + rowSums(res2$lambda_component[1:5, c("lambda11", "lambda12")]) res2$lambda[1:5, "lambda1"] ## ----------------------------------------------------------------------------- inter_arrival2 <- res2$inter_arrival type2 <- res2$type ## ----------------------------------------------------------------------------- logLik(hspec2, inter_arrival = inter_arrival2, type = type2) ## ---- warning=FALSE----------------------------------------------------------- hspec0 <- hspec2 mle <- hfit(hspec0, inter_arrival = inter_arrival2, type = type2) summary(mle) ## ----------------------------------------------------------------------------- coef(mle) ## ----------------------------------------------------------------------------- miscTools::stdEr(mle) ## ----------------------------------------------------------------------------- print(hspec2) ## ---- warning=FALSE----------------------------------------------------------- res2 <- hsim(hspec2, size = 1000) ## ---- warning=FALSE----------------------------------------------------------- mu0 <- matrix(c(0.15, 0.15), nrow = 2) alpha0 <- matrix(c(0.75, 0.75, 0.75, 0.75), nrow = 2, byrow=TRUE) beta0 <- matrix(c(2.6, 2.6, 2.6, 2.6), nrow = 2, byrow=TRUE) hspec0 <- new("hspec", mu=mu0, alpha=alpha0, beta=beta0) summary(hfit(hspec0, inter_arrival = res2$inter_arrival, type = res2$type)) ## ---- warning=FALSE----------------------------------------------------------- mu0 <- matrix(c(0.15, 0.15), nrow = 2) alpha0 <- matrix(c(0.75, 0.751, 0.751, 0.75), nrow = 2, byrow=TRUE) beta0 <- matrix(c(2.6, 2.6, 2.6, 2.6), nrow = 2, byrow=TRUE) hspec0 <- new("hspec", mu=mu0, alpha=alpha0, beta=beta0) summary(hfit(hspec0, inter_arrival = res2$inter_arrival, type = res2$type)) ## ---- warning=FALSE----------------------------------------------------------- mu0 <- matrix(c(0.15, 0.14), nrow = 2) alpha0 <- matrix(c(0.75, 0.751, 0.751, 0.75), nrow = 2, byrow=TRUE) beta0 <- matrix(c(2.6, 2.6, 2.6, 2.6), nrow = 2, byrow=TRUE) hspec0 <- new("hspec", mu=mu0, alpha=alpha0, beta=beta0) summary(hfit(hspec0, inter_arrival = res2$inter_arrival, type = res2$type)) ## ---- warning=FALSE----------------------------------------------------------- summary(hfit(hspec2, inter_arrival = res2$inter_arrival, type = res2$type, reduced=FALSE)) ## ---- warning=FALSE, message = FALSE, fig.height=5, fig.width=5--------------- hrp <- new("hspec", mu = 0.3, alpha = 1.2, beta = 1.5) res_rp <- hsim(hrp, size = 1000) rp <- residual_process(component = 1, inter_arrival = res_rp$inter_arrival, type = res_rp$type, rambda_component = res_rp$rambda_component, mu = 0.3, beta = 1.5) p <- ppoints(100) q <- quantile(rp,p=p) plot(qexp(p), q, xlab="Theoretical Quantiles",ylab="Sample Quantiles") qqline(q, distribution=qexp,col="blue", lty=2) ## ---- message=FALSE, fig.height=5, fig.width=5-------------------------------- # estimation mle_rp <- hfit(new("hspec", mu = 0.2, alpha = 1, beta = 2), res_rp$inter_arrival) # construct hspec from estimation result he <- new("hspec", mu = coef(mle_rp)["mu1"], alpha = coef(mle_rp)["alpha1"], beta = coef(mle_rp)["beta1"]) # infer intensity infered_res <- infer_lambda(he, res_rp$inter_arrival, res_rp$type) # compute residuals where we use rp2 <- residual_process(component = 1, inter_arrival = res_rp$inter_arrival, type = res_rp$type, rambda_component = infered_res$rambda_component, mu = coef(mle_rp)["mu1"], beta = coef(mle_rp)["beta1"]) p <- ppoints(100) q <- quantile(rp2, p=p) plot(qexp(p), q, xlab="Theoretical Quantiles",ylab="Sample Quantiles") qqline(q, distribution=qexp,col="blue", lty=2) ## ----------------------------------------------------------------------------- mu <- matrix(c(0.02, 0.02), nrow=2) beta_1 <- matrix(rep(10, 4), nrow=2) beta_2 <- matrix(rep(1, 4), nrow=2) beta <- cbind(beta_1, beta_2) alpha_1 <- matrix(c(3, 2, 2, 3), nrow=2, byrow=TRUE) alpha_2 <- matrix(c(0.3, 0.2, 0.2, 0.3), nrow=2, byrow=TRUE) alpha <- cbind(alpha_1, alpha_2) print(alpha) ## ----------------------------------------------------------------------------- type_col_map <- list(c(1,3), # columns for dN1 c(2,4)) # columns for dN2 type_col_map ## ----------------------------------------------------------------------------- cat("Part of alpha associated with N1: \n") alpha[, type_col_map[[1]]] # associated with N1 cat("Part of alpha associated with N2: \n") alpha[, type_col_map[[2]]] # associated with N2 cat("Part of beta associated with N1: \n") beta[, type_col_map[[1]]] # associated with N1 cat("Part of beta associated with N2: \n") beta[, type_col_map[[2]]] # associated with N2 ## ----------------------------------------------------------------------------- h <- new("hspec", mu = mu, alpha = alpha, beta=beta, type_col_map = type_col_map) h ## ----------------------------------------------------------------------------- res_mk <- hsim(h, size = 2000, # for an illustration purpose lambda_component0 = matrix(seq(1, 1.7, 0.1), nrow = 2)) res_mk ## ---- warning=FALSE----------------------------------------------------------- summary(hfit(h, res_mk$inter_arrival, res_mk$type, lambda_component0 = matrix(seq(1, 1.7, 0.1), nrow = 2))) ## ----------------------------------------------------------------------------- mu <- function(param = c(theta_p = 0.15, theta_n = 0.21, kappa = 0.12)){ theta <- (param["theta_n"] + param["theta_p"])/2 theta_tl <- (param["theta_n"] - param["theta_p"])/2 matrix(c(theta/2/(1 - param["kappa"]) + theta_tl/2/(1 + param["kappa"]), theta/2/(1 - param["kappa"]) - theta_tl/2/(1 + param["kappa"])), nrow=2) } alpha <- function(param = c(eta = 5, nu = 3)){ zeta <- (param["eta"] + param["nu"])/2 zeta_tl <- (param["eta"] - param["nu"])/2 matrix(c(zeta, zeta_tl, zeta, -zeta_tl, zeta, -zeta_tl, zeta, zeta_tl), nrow=2, byrow=TRUE) } beta <- function(param = c(beta = 12, kappa = 0.12)){ beta1 <- param["beta"] * (1 - param["kappa"]) beta2 <- param["beta"] * (1 + param["kappa"]) matrix(c(beta1, beta2, beta1, beta2, beta1, beta2, beta1, beta2), nrow = 2, byrow = TRUE) } type_col_map <- list(c(1,2), c(3,4)) h_sy <- new("hspec", mu = mu, alpha = alpha, beta = beta, type_col_map = type_col_map) h_sy ## ----------------------------------------------------------------------------- # run simulation res_sy <- hsim(h_sy, size = 2000, lambda_component0 = matrix(rep(1, 2 * 4), nrow=2)) summary(res_sy) ## ---- warning=FALSE----------------------------------------------------------- fit_sy <- hfit(h_sy, inter_arrival=res_sy$inter_arrival, type=res_sy$type, lambda_component0 = matrix(rep(1, 2 * 4), nrow=2)) summary(fit_sy) ## ---- warning=FALSE----------------------------------------------------------- mu <- matrix(c(0.15, 0.15), nrow=2) alpha <- matrix(c(0.75, 0.6, 0.6, 0.75), nrow=2, byrow=T) beta <- matrix(c(2.6, 2.6, 2.6, 2.6), nrow=2) rmark <- function(param = c(p=0.65), ...){ rgeom(1, p=param[1]) + 1 } impact <- function(param = c(eta1=0.2), alpha, n, mark, ...){ ma <- matrix(rep(mark[n]-1, 4), nrow = 2) ma * matrix( rep(param["eta1"], 4), nrow=2) } hi <- new("hspec", mu=mu, alpha=alpha, beta=beta, rmark = rmark, impact=impact) hi ## ---- warning=FALSE----------------------------------------------------------- res_impact <- hsim(hi, size=1000, lambda_component0 = matrix(rep(0.1,4), nrow=2)) summary(res_impact) ## ---- warning=FALSE----------------------------------------------------------- fit <- hfit(hi, inter_arrival = res_impact$inter_arrival, type = res_impact$type, mark = res_impact$mark, lambda_component0 = matrix(rep(0.1,4), nrow=2)) summary(fit) ## ----------------------------------------------------------------------------- rmark <- function(param = c(p=0.65), ...){ rgeom(1, p=param[1]) + 1 } h <- new("hspec", mu=0.15, alpha=0.7, beta=1.6, eta=0.3, rmark = rmark) h ## ----------------------------------------------------------------------------- res <- hsim(h, size = 1000) summary(res) ## ---- warning=FALSE----------------------------------------------------------- fit <- hfit(h, inter_arrival = res$inter_arrival, type = res$type, mark = res$mark) summary(fit) ## ---- warning=FALSE----------------------------------------------------------- h_md <- h h_md@dmark <- function(param = c(p = 0.1), n=n, mark=mark, ...){ dgeom(mark[n] - 1, prob = param["p"]) } mle_md <- hfit(h_md, inter_arrival = res$inter_arrival, type = res$type, mark = res$mark) summary(mle_md) ## ---- warning=FALSE----------------------------------------------------------- mu <- matrix(c(0.02, 0.02, 0.04, 0.04), nrow = 4) alpha <- function(param = c(alpha11 = 0.2, alpha12 = 0.3, alpha33 = 0.3, alpha34 = 0.4)){ matrix(c(param["alpha11"], param["alpha12"], 0, 0, param["alpha12"], param["alpha11"], 0, 0, 0, 0, param["alpha33"], param["alpha34"], 0, 0, param["alpha34"], param["alpha33"]), nrow = 4, byrow = TRUE) } beta <- matrix(c(rep(0.7, 8), rep(1.1, 8)), nrow = 4, byrow = TRUE) ## ----------------------------------------------------------------------------- impact <- function(param = c(alpha1n=0.25, alpha1w=0.1, alpha2n=0.1, alpha2w=0.2), N=N, n=n, ...){ Psi <- matrix(c(0, 0, param['alpha1w'], param['alpha1n'], 0, 0, param['alpha1n'], param['alpha1w'], param['alpha2w'], param['alpha2n'], 0, 0, param['alpha2n'], param['alpha2w'], 0, 0), nrow=4, byrow=TRUE) ind <- N[,"N1"][n] - N[,"N2"][n] > N[,"N3"][n] - N[,"N4"][n] km <- matrix(c(!ind, !ind, !ind, !ind, ind, ind, ind, ind, ind, ind, ind, ind, !ind, !ind, !ind, !ind), nrow = 4, byrow = TRUE) km * Psi } hspec_fl <- new("hspec", mu = mu, alpha = alpha, beta = beta, impact = impact) hspec_fl ## ----------------------------------------------------------------------------- hr_fl <- hsim(hspec_fl, size=1000) summary(hr_fl) ## ---- warning=FALSE----------------------------------------------------------- fit_fl <- hfit(hspec_fl, hr_fl$inter_arrival, hr_fl$type) summary(fit_fl) ## ----------------------------------------------------------------------------- # presumed initial bid and ask prices initial_ask_price <- 1250 #cents initial_bid_price <- 1150 #cents initial_level <- round((initial_ask_price - initial_bid_price) - 1) initial_mid_price <- (initial_bid_price + initial_ask_price) / 2 mu <- function(param = c(mu1 = 0.08, zeta1 = 0.10), n=n, Nc=Nc, ...){ if(n == 1){ level <- initial_level mid <- initial_mid_price } else { level <- Nc[n-1,1] - Nc[n-1,2] - (Nc[n-1,3] - Nc[n-1,4]) + initial_level ask <- initial_ask_price + (Nc[n-1,1] - Nc[n-1,2]) bid <- initial_bid_price + (Nc[n-1,3] - Nc[n-1,4]) mid <- (ask + bid) / 2 } if(level <= 0){ matrix(c(param["mu1"], 0, 0, param["mu1"]), nrow = 4) } else { matrix(c(param["mu1"], param["zeta1"] * level / mid, param["zeta1"]*level / mid, param["mu1"]), nrow = 4) } } ## ---- warning=FALSE----------------------------------------------------------- alpha <- function(param = c(alpha_s1=4, alpha_m=26, alpha_s2=5, alpha_w1=11, alpha_w2=7)){ matrix(c(param["alpha_s1"], param["alpha_m"], param["alpha_s2"], 0, param["alpha_w1"], 0, 0, param["alpha_w2"], param["alpha_w2"], 0, 0, param["alpha_w1"], 0, param["alpha_s2"], param["alpha_m"], param["alpha_s1"]), nrow = 4, byrow = TRUE) } impact <- function(param = c(xi = 2.7), n=n, Nc=Nc, lambda_component = lambda_component, ... ){ if(n == 1){ level <- initial_level # mid <- initial_mid_price } else { level <- Nc[n,1] - Nc[n,2] - (Nc[n,3] - Nc[n,4]) + initial_level ask <- initial_ask_price + (Nc[n,1] - Nc[n,2]) bid <- initial_bid_price + (Nc[n,3] - Nc[n,4]) mid <- (ask + bid)/2 } lambda_component_matrix <- matrix(lambda_component[n, ], nrow=4, byrow=TRUE) l2 <- sum(lambda_component_matrix[2,]) # sum of second row l3 <- sum(lambda_component_matrix[3,]) # sum of third row im <- matrix(c(0, 0, 0, 0, 0, -l2 + param["xi"]*level/mid, -l2 + param["xi"]*level/mid, 0, 0, -l3 + param["xi"]*level/mid, -l3 + param["xi"]*level/mid, 0, 0, 0, 0, 0), nrow = 4, byrow = TRUE) } beta <- matrix(rep(50, 16), nrow = 4, byrow=TRUE) rmark <- function(n=n, Nc=Nc, type, ...){ if(n == 1){ level <- initial_level } else { level <- Nc[n-1,1] - Nc[n-1,2] - (Nc[n-1,3] - Nc[n-1,4]) + initial_level } if (type[n] == 2 | type[n] == 3){ min(level, rgeom(1, p=0.65) + 1) } else { rgeom(1, p=0.65) + 1 } } h_ba <- new("hspec", mu = mu, alpha = alpha, beta = beta, impact=impact, rmark = rmark) h_ba ## ----------------------------------------------------------------------------- hr_ba <- hsim(h_ba, size=1000, lambda_component0 = matrix(rep(1, 16), 4)) summary(hr_ba) ## ---- warning=FALSE----------------------------------------------------------- mle_ba <- hfit(h_ba, inter_arrival = hr_ba$inter_arrival, type = hr_ba$type, lambda_component0 = matrix(rep(1, 16), 4)) summary(mle_ba)