1efp <- function(formula, data = list(), 2 type = c("Rec-CUSUM", "OLS-CUSUM", "Rec-MOSUM", "OLS-MOSUM", 3 "RE", "ME", "Score-CUSUM", "Score-MOSUM", "fluctuation"), 4 h = 0.15, dynamic = FALSE, rescale = TRUE, lrvar = FALSE, vcov = NULL) 5{ 6 if(!inherits(formula, "formula")) { 7 mt <- terms(formula) 8 X <- if(is.matrix(formula$x)) 9 formula$x 10 else model.matrix(mt, model.frame(formula)) 11 y <- if(is.vector(formula$y)) 12 formula$y 13 else model.response(model.frame(formula)) 14 } else { 15 mf <- model.frame(formula, data = data) 16 mt <- attr(mf, "terms") 17 y <- model.response(mf) 18 X <- model.matrix(formula, data = data) 19 } 20 21 lmfit <- function(x, y, ...) { 22 rval <- lm.fit(x, y, ...) 23 rval$terms <- mt 24 rval$x <- x 25 class(rval) <- "lm" 26 return(rval) 27 } 28 29 lrvartype <- if(is.character(lrvar)) lrvar else "Andrews" 30 sdev <- if(identical(lrvar, FALSE)) { 31 function(x, df = NULL) { 32 if(is.null(df)) df <- length(x) - 1 33 sd(x) * sqrt((length(x) - 1)/df) 34 } 35 } else { 36 function(x, df = NULL) { 37 s <- sqrt(sandwich::lrvar(x, type = lrvartype) * length(x)) 38 if(!is.null(df)) s <- s * sqrt(length(x)/df) 39 return(s) 40 } 41 } 42 43 n <- nrow(X) 44 if(dynamic) 45 { 46 Xnames <- colnames(X) 47 X <- cbind(y[1:(n-1)],X[2:n,]) 48 colnames(X) <- c("lag", Xnames) 49 y <- y[-1] 50 n <- n-1 51 } 52 k <- ncol(X) 53 type <- match.arg(type) 54 if(type == "fluctuation") type <- "RE" 55 56 retval <- list(process = NULL, 57 type = type, 58 nreg = k, 59 nobs = n, 60 call = match.call(), 61 formula = formula, 62 par = NULL, 63 type.name = NULL, 64 lim.process = NULL, 65 coef = NULL, 66 Q12 = NULL, 67 datatsp = NULL, 68 rescale = rescale) 69 70 orig.y <- NULL 71 72 switch(type, 73 74 ## empirical process of Recursive CUSUM model 75 76 "Rec-CUSUM" = { 77 w <- recresid(X, y) 78 sigma <- sdev(w) ## sqrt(var(w)) 79 process <- cumsum(c(0,w))/(sigma*sqrt(n-k)) 80 if(is.ts(data)) { 81 if(NROW(data) == n) process <- ts(process, end = end(data), frequency = frequency(data)) 82 } else { 83 env <- environment(formula) 84 if(missing(data)) data <- env 85 orig.y <- eval(attr(mt, "variables")[[2]], data, env) 86 if(is.ts(orig.y) && (NROW(orig.y) == n)) 87 process <- ts(process, end = end(orig.y), 88 frequency = frequency(orig.y)) 89 } 90 retval$type.name <- "Recursive CUSUM test" 91 retval$lim.process <- "Brownian motion" 92 }, 93 94 ## empirical process of OLS-based CUSUM model 95 96 "OLS-CUSUM" = { 97 fm <- lm.fit(X,y) 98 e <- fm$residuals 99 sigma <- sdev(e, df = fm$df.residual) ## sqrt(sum(e^2)/fm$df.residual) 100 process <- cumsum(c(0,e))/(sigma * sqrt(n)) 101 if(is.ts(data)) { 102 if(NROW(data) == n) process <- ts(process, end = end(data), frequency = frequency(data)) 103 } else { 104 env <- environment(formula) 105 if(missing(data)) data <- env 106 orig.y <- eval(attr(mt, "variables")[[2]], data, env) 107 if(is.ts(orig.y) && (NROW(orig.y) == n)) 108 process <- ts(process, end = end(orig.y), 109 frequency = frequency(orig.y)) 110 } 111 retval$type.name <- "OLS-based CUSUM test" 112 retval$lim.process <- "Brownian bridge" 113 }, 114 115 ## empirical process of Recursive MOSUM model 116 117 "Rec-MOSUM" = { 118 w <- recresid(X, y) 119 nw <- n - k 120 nh <- floor(nw*h) 121 process <- rep(0, (nw-nh)) 122 for(i in 0:(nw-nh)) 123 { 124 process[i+1] <- sum(w[(i+1):(i+nh)]) 125 } 126 sigma <- sdev(w, df = nw - k) ## sqrt(var(w)*(nw-1)/(nw-k)) 127 process <- process/(sigma * sqrt(nw)) 128 if(is.ts(data)) { 129 if(NROW(data) == n) process <- ts(process, end = time(data)[(n-floor(0.5 + nh/2))], frequency = frequency(data)) 130 } else { 131 env <- environment(formula) 132 if(missing(data)) data <- env 133 orig.y <- eval(attr(mt, "variables")[[2]], data, env) 134 if(is.ts(orig.y) && (NROW(orig.y) == n)) { 135 process <- ts(process, end = time(orig.y)[(n-floor(0.5 + nh/2))], 136 frequency = frequency(orig.y)) 137 } else { 138 process <- ts(process, end = (n-floor(0.5 + nh/2))/n, 139 frequency = n) 140 } 141 } 142 retval$par <- h 143 retval$type.name <- "Recursive MOSUM test" 144 retval$lim.process <- "Brownian motion increments" 145 }, 146 147 ## empirical process of OLS-based MOSUM model 148 149 "OLS-MOSUM" = { 150 fm <- lm.fit(X,y) 151 e <- fm$residuals 152 sigma <- sdev(e, df = fm$df.residual) ## sqrt(sum(e^2)/fm$df.residual) 153 nh <- floor(n*h) 154 process <- cumsum(c(0,e)) 155 process <- process[-(1:nh)] - process[1:(n-nh+1)] 156 process <- process/(sigma * sqrt(n)) 157 if(is.ts(data)) { 158 if(NROW(data) == n) process <- ts(process, end = time(data)[(n-floor(0.5 + nh/2))], frequency = frequency(data)) 159 } else { 160 env <- environment(formula) 161 if(missing(data)) data <- env 162 orig.y <- eval(attr(mt, "variables")[[2]], data, env) 163 if(is.ts(orig.y) && (NROW(orig.y) == n)) { 164 process <- ts(process, end = time(orig.y)[(n-floor(0.5 + nh/2))], 165 frequency = frequency(orig.y)) 166 } else { 167 process <- ts(process, end = (n-floor(0.5 + nh/2))/n, 168 frequency = n) 169 } 170 } 171 retval$par <- h 172 retval$type.name <- "OLS-based MOSUM test" 173 retval$lim.process <- "Brownian bridge increments" 174 }, 175 176 ## empirical process of recursive estimates fluctuation 177 178 "RE" = { 179 m.fit <- lmfit(X,y) 180 beta.hat <- m.fit$coefficients 181 sigma <- sdev(m.fit$residual, df = m.fit$df.residual) ## sqrt(sum(m.fit$residual^2)/m.fit$df.residual) 182 process <- matrix(rep(0,k*(n-k+1)), nrow=k) 183 Q12 <- if(is.null(vcov)) { 184 root.matrix(crossprod(X))/(sigma * sqrt(n)) 185 } else { 186 root.matrix(solve(vcov(m.fit)))/sqrt(n) 187 } 188 if(rescale) { 189 for(i in k:(n-1)) { 190 mi.fit <- lmfit(as.matrix(X[1:i,]), y[1:i]) 191 if(!is.null(vcov)) warning("custom vcov not implemented yet if rescale = TRUE") 192 Qi12 <- root.matrix(crossprod(X[1:i, ]))/(sigma * sqrt(i)) 193 process[,(i-k+1)] <- Qi12 %*% (mi.fit$coefficients - beta.hat) 194 } 195 } else { 196 for(i in k:(n-1)) { 197 mi.fit <- lmfit(as.matrix(X[1:i,]), y[1:i]) 198 process[,(i-k+1)] <- Q12 %*% (mi.fit$coefficients - beta.hat) 199 } 200 } 201 process <- t(cbind(0, process)) * matrix(rep((k - 1):n, k), ncol = k)/sqrt(n) 202 colnames(process) <- colnames(X) 203 if(is.ts(data)) { 204 if(NROW(data) == n) process <- ts(process, end = end(data), frequency = frequency(data)) 205 } else { 206 env <- environment(formula) 207 if(missing(data)) data <- env 208 orig.y <- eval(attr(mt, "variables")[[2]], data, env) 209 if(is.ts(orig.y) && (NROW(orig.y) == n)) { 210 process <- ts(process, end = end(orig.y), 211 frequency = frequency(orig.y)) 212 } 213 } 214 retval$Q12 <- Q12 215 retval$type.name <- "RE test (recursive estimates test)" 216 retval$lim.process <- "Brownian bridge" 217 }, 218 219 ## empirical process of moving estimates fluctuation 220 221 "ME" = { 222 m.fit <- lmfit(X,y) 223 beta.hat <- m.fit$coefficients 224 sigma <- sdev(m.fit$residual, df = m.fit$df.residual) ## sqrt(sum(m.fit$residual^2)/m.fit$df.residual) 225 nh <- floor(n*h) 226 process <- matrix(rep(0,k*(n-nh+1)), nrow=k) 227 Q12 <- if(is.null(vcov)) { 228 root.matrix(crossprod(X))/(sigma * sqrt(n)) 229 } else { 230 root.matrix(solve(vcov(m.fit)))/sqrt(n) 231 } 232 if(rescale) { 233 for(i in 0:(n-nh)) { 234 mi.fit <- lmfit(as.matrix(X[(i+1):(i+nh),]), y[(i+1):(i+nh)]) 235 Qnh12 <- if(is.null(vcov)) { 236 root.matrix(crossprod(X[(i+1):(i+nh),]))/(sigma * sqrt(nh)) 237 } else { 238 root.matrix(solve(vcov(mi.fit)))/sqrt(nh) 239 } 240 process[, i+1] <- Qnh12 %*% (mi.fit$coefficients - beta.hat) 241 } 242 } else { 243 for(i in 0:(n-nh)) { 244 mi.fit <- lmfit(as.matrix(X[(i+1):(i+nh),]), y[(i+1):(i+nh)]) 245 process[, i+1] <- Q12 %*% (mi.fit$coefficients - beta.hat) 246 } 247 } 248 process <- nh * t(process)/sqrt(n) 249 colnames(process) <- colnames(X) 250 if(is.ts(data)) { 251 if(NROW(data) == n) process <- ts(process, end = time(data)[(n-floor(0.5 + nh/2))], frequency = frequency(data)) 252 } else { 253 env <- environment(formula) 254 if(missing(data)) data <- env 255 orig.y <- eval(attr(mt, "variables")[[2]], data, env) 256 if(is.ts(orig.y) && (NROW(orig.y) == n)) { 257 process <- ts(process, end = time(orig.y)[(n-floor(0.5 + nh/2))], 258 frequency = frequency(orig.y)) 259 } else { 260 process <- ts(process, end = (n-floor(0.5 + nh/2))/n, 261 frequency = n) 262 } 263 } 264 retval$par <- h 265 retval$Q12 <- Q12 266 retval$type.name <- "ME test (moving estimates test)" 267 retval$lim.process <- "Brownian bridge increments" 268 }, 269 270 "Score-CUSUM" = { 271 fm <- lm.fit(X,y) 272 e <- as.vector(fm$residuals) 273 sigma2 <- sum(e^2)/n 274 k <- k + 1 275 ## Q12 <- sqrt(sigma2) * root.matrix(crossprod(X))/sqrt(n) 276 ## Q12 <- rbind(cbind(Q12, 0), 0) 277 ## Q12[k,k] <- sqrt(2)*sigma2 278 279 process <- cbind(X * e, e^2 - sigma2)/sqrt(n) 280 Q12 <- root.matrix(crossprod(process)) 281 process <- rbind(0, process) 282 process <- apply(process, 2, cumsum) 283 process <- t(chol2inv(chol(Q12)) %*% t(process)) 284 285 colnames(process) <- c(names(coef(fm)), "(Variance)") 286 if(is.ts(data)) { 287 if(NROW(data) == n) process <- ts(process, end = end(data), frequency = frequency(data)) 288 } else { 289 env <- environment(formula) 290 if(missing(data)) data <- env 291 orig.y <- eval(attr(mt, "variables")[[2]], data, env) 292 if(is.ts(orig.y) && (NROW(orig.y) == n)) { 293 process <- ts(process, end = end(orig.y), 294 frequency = frequency(orig.y)) 295 } 296 } 297 retval$type.name <- "Score-based CUSUM test" 298 retval$lim.process <- "Brownian bridge" 299 retval$Q12 <- Q12 300 }, 301 302 "Score-MOSUM" = { 303 fm <- lm.fit(X,y) 304 e <- as.vector(fm$residuals) 305 sigma2 <- sum(e^2)/n 306 k <- k + 1 307 nh <- floor(n*h) 308 ## Q12 <- sqrt(sigma2) * root.matrix(crossprod(X))/sqrt(n) 309 ## Q12 <- rbind(cbind(Q12, 0), 0) 310 ## Q12[k,k] <- sqrt(2)*sigma2 311 312 process <- cbind(X * e, e^2 - sigma2)/sqrt(n) 313 Q12 <- root.matrix(crossprod(process)) 314 process <- rbind(0, process) 315 process <- apply(process, 2, cumsum) 316 process <- process[-(1:nh),] - process[1:(n-nh+1),] 317 process <- t(chol2inv(chol(Q12)) %*% t(process)) 318 319 colnames(process) <- c(names(coef(fm)), "(Variance)") 320 if(is.ts(data)) { 321 if(NROW(data) == n) process <- ts(process, end = time(data)[(n-floor(0.5 + nh/2))], frequency = frequency(data)) 322 } else { 323 env <- environment(formula) 324 if(missing(data)) data <- env 325 orig.y <- eval(attr(mt, "variables")[[2]], data, env) 326 if(is.ts(orig.y) && (NROW(orig.y) == n)) { 327 process <- ts(process, end = time(orig.y)[(n-floor(0.5 + nh/2))], 328 frequency = frequency(orig.y)) 329 } else { 330 process <- ts(process, end = (n-floor(0.5 + nh/2))/n, 331 frequency = n) 332 } 333 } 334 retval$par <- h 335 retval$type.name <- "Score-based MOSUM test" 336 retval$lim.process <- "Brownian bridge increments" 337 retval$Q12 <- Q12 338 }) 339 340 341 if(!is.ts(process)) 342 process <- ts(process, start = 0, frequency = (NROW(process)-1)) 343 344 retval$process <- process 345 346 if(is.ts(data) & NROW(data) == n) 347 retval$datatsp <- tsp(data) 348 else if(!is.null(orig.y) && is.ts(orig.y) & NROW(orig.y) == n) 349 retval$datatsp <- tsp(orig.y) 350 else 351 retval$datatsp <- c(0, 1, n) 352 353 m.fit <- lm.fit(X,y) 354 retval$coefficients <- coefficients(m.fit) 355 retval$sigma <- sqrt(sum(m.fit$residual^2)/m.fit$df.residual) 356 class(retval) <- c("efp") 357 return(retval) 358} 359 360 361plot.efp <- function(x, alpha = 0.05, alt.boundary = FALSE, boundary = TRUE, 362 functional = "max", main = NULL, ylim = NULL, 363 ylab = "Empirical fluctuation process", ...) 364{ 365 if(is.null(functional)) fun <- "max" 366 else fun <- match.arg(functional, c("max", "range", "maxL2", "meanL2")) 367 bound <- boundary(x, alpha = alpha, alt.boundary = alt.boundary, functional = fun) 368 pos <- FALSE 369 ave <- FALSE 370 371 if(is.null(main)){ 372 if(alt.boundary & fun == "max" & (x$lim.process %in% c("Brownian motion", "Brownian bridge"))){ 373 main <- paste(x$type.name, "with alternative boundaries") 374 } 375 else { 376 if(alt.boundary) warning("no alternative boundaries available") 377 if(fun == "meanL2") 378 main <- paste(x$type.name, "with mean L2 norm") 379 else if(fun == "maxL2") 380 main <- paste(x$type.name, "with max L2 norm") 381 else 382 main <- x$type.name 383 } 384 } 385 386 if(!is.null(functional) && x$lim.process %in% c("Brownian bridge", "Brownian bridge increments")) { 387 z <- as.matrix(x$process) 388 k <- ncol(z) 389 390 switch(functional, 391 "max" = { 392 if(k > 1) { 393 z <- apply(abs(z), 1, max) 394 pos <- TRUE 395 } 396 }, 397 "range" = { stop("no plot available for range functional") }, 398 "maxL2" = { 399 if(x$lim.process == "Brownian bridge") { 400 z <- rowSums(z^2) 401 pos <- TRUE 402 } else { 403 stop("no test/plot available for mean L2 functional") 404 } 405 }, 406 407 "meanL2" = { 408 if(x$lim.process == "Brownian bridge") { 409 z <- rowSums(z^2) 410 ave <- TRUE 411 pos <- TRUE 412 } else { 413 stop("no test/plot available for mean L2 functional") 414 } 415 }) 416 z <- ts(as.vector(z), start = start(x$process), frequency = frequency(x$process)) 417 } else { 418 z <- x$process 419 } 420 421 if(is.null(ylim)) { 422 ymax <- max(c(z, bound)) 423 if(pos) ymin <- 0 424 else ymin <- min(c(z, -bound)) 425 ylim <- c(ymin, ymax) 426 } 427 428 if(boundary) 429 panel <- function(y, ...) { 430 lines(y, ...) 431 lines(bound, col=2) 432 lines(-bound, col=2) 433 abline(0,0) 434 } 435 else 436 panel <- function(y, ...) { 437 lines(y, ...) 438 abline(0,0) 439 } 440 441 if(any(attr(z, "class") == "mts")) 442 plot(z, main = main, panel = panel, ...) 443 else { 444 plot(z, main = main, ylab = ylab, ylim = ylim, ...) 445 if(boundary) { 446 lines(bound, col=2) 447 if(!pos) lines(-bound, col=2) 448 if(ave) { 449 avez <- ts(rep(mean(z), length(bound)), start = start(bound), frequency = frequency(bound)) 450 lines(avez, lty = 2) 451 } 452 } 453 abline(0,0) 454 } 455} 456 457pvalue.efp <- function(x, lim.process, alt.boundary, functional = "max", h = NULL, k = NULL) 458{ 459 lim.process <- match.arg(lim.process, 460 c("Brownian motion", "Brownian bridge", "Brownian motion increments", "Brownian bridge increments")) 461 functional <- match.arg(functional, c("max", "range", "meanL2", "maxL2")) 462 463 switch(lim.process, 464 465 "Brownian motion" = { 466 if(functional == "max") { 467 if(alt.boundary) 468 { 469 pval <- c(1, 0.997, 0.99, 0.975, 0.949, 0.912, 0.864, 0.806, 0.739, 0.666, 470 0.589, 0.512, 0.437, 0.368, 0.307, 0.253, 0.205, 0.163, 0.129, 0.100, 471 0.077, 0.058, 0.043, 0.032, 0.024, 0.018, 0.012, 0.009, 0.006, 0.004, 472 0.003, 0.002, 0.001, 0.001, 0.001) 473 critval <- (10:44)/10 474 p <- approx(critval, pval, x, rule=2)$y 475 } else { 476 p <- ifelse(x < 0.3, 1 - 0.1465*x, 477 2*(1-pnorm(3*x) + exp(-4*(x^2))*(pnorm(x)+pnorm(5*x)-1)-exp(-16*(x^2))*(1-pnorm(x)))) 478 } 479 } else { 480 stop("only max functional implemented for Brownian motion") 481 } 482 p <- 1 - (1-p)^k 483 }, 484 485 "Brownian bridge" = { 486 switch(functional, 487 488 "max" = { 489 if(alt.boundary) 490 { 491 pval <- c(1, 1, 0.997, 0.99, 0.977, 0.954, 0.919, 0.871, 0.812, 0.743, 492 0.666, 0.585, 0.504, 0.426, 0.353, 0.288, 0.230, 0.182, 0.142, 0.109, 493 0.082, 0.062, 0.046, 0.034, 0.025, 0.017, 0.011, 0.008, 0.005, 0.004, 494 0.003, 0.002, 0.001, 0.001, 0.0001) 495 critval <- (12:46)/10 496 p <- approx(critval, pval, x, rule=2)$y 497 p <- 1 - (1-p)^k 498 } else { 499 p <- ifelse(x<0.1, 1, 500 { 501 summand <- function(a,b) 502 { 503 exp(-2*(a^2)*(b^2))*(-1)^a 504 } 505 p <- 0 506 for(i in 1:100) p <- p + summand(i,x) 507 1-(1+2*p)^k 508 }) 509 } 510 }, 511 512 "range" = { 513 p <- ifelse(x<0.4,1, 514 { 515 p <- 0 516 for(i in 1:10) p <- p + (4*i^2*x^2 - 1) * exp(-2*i^2*x^2) 517 1-(1-2*p)^k 518 }) 519 }, 520 521 "maxL2" = { 522 if(k > 25) { 523 warning("number of regressors > 25, critical values for 25 regressors used") 524 k <- 25 525 } 526 critval <- get("sc.maxL2")[as.character(k), ] 527 p <- approx(c(0, critval), c(1, 1-as.numeric(names(critval))), x, rule=2)$y 528 }, 529 530 "meanL2" = { 531 if(k > 25) { 532 warning("number of regressors > 25, critical values for 25 regressors used") 533 k <- 25 534 } 535 critval <- get("sc.meanL2")[as.character(k), ] 536 p <- approx(c(0, critval), c(1, 1-as.numeric(names(critval))), x, rule=2)$y 537 }) 538 }, 539 540 "Brownian motion increments" = { 541 if(alt.boundary) stop("alternative boundaries for Brownian motion increments not available") 542 if(functional != "max") stop("only max functional for Brownian motion increments available") 543 544 crit.table <- cbind(c(3.2165, 2.9795, 2.8289, 2.7099, 2.6061, 2.5111, 2.4283, 2.3464, 2.2686, 2.2255), 545 c(3.3185, 3.0894, 2.9479, 2.8303, 2.7325, 2.6418, 2.5609, 2.4840, 2.4083, 2.3668), 546 c(3.4554, 3.2368, 3.1028, 2.9874, 2.8985, 2.8134, 2.7327, 2.6605, 2.5899, 2.5505), 547 c(3.6622, 3.4681, 3.3382, 3.2351, 3.1531, 3.0728, 3.0043, 2.9333, 2.8743, 2.8334), 548 c(3.8632, 3.6707, 3.5598, 3.4604, 3.3845, 3.3102, 3.2461, 3.1823, 3.1229, 3.0737), 549 c(4.1009, 3.9397, 3.8143, 3.7337, 3.6626, 3.5907, 3.5333, 3.4895, 3.4123, 3.3912)) 550 tablen <- dim(crit.table)[2] 551 tableh <- (1:10)*0.05 552 tablep <- c(0.2, 0.15, 0.1, 0.05, 0.025, 0.01) 553 554 ## crit.table gives Table 1 of Chu, Hornik, Kuan (1995) 555 ## but the corresponding test statistic is scaled differently 556 ## by a factor of sqrt(h). 557 crit.table <- crit.table * sqrt(tableh) 558 559 tableipl <- numeric(tablen) 560 for(i in (1:tablen)) tableipl[i] <- approx(tableh, crit.table[,i], h, rule = 2)$y 561 p <- approx(c(0,tableipl), c(1,tablep), x, rule = 2)$y 562 p <- 1 - (1-p)^k 563 }, 564 565 "Brownian bridge increments" = { 566 if(k > 6) k <- 6 567 switch(functional, 568 "max" = { 569 crit.table <- get("sc.me")[(((k-1)*10+1):(k*10)), ] 570 tablen <- dim(crit.table)[2] 571 tableh <- (1:10)*0.05 572 tablep <- c(0.1, 0.05, 0.025, 0.01) 573 tableipl <- numeric(tablen) 574 for(i in (1:tablen)) tableipl[i] <- approx(tableh, crit.table[,i], h, rule = 2)$y 575 p <- approx(c(0,tableipl), c(1,tablep), x, rule = 2)$y 576 }, 577 "range" = { 578 p <- ifelse(x<0.53,1, 579 { 580 p <- 0 581 for(i in 1:10) p <- p + (-1)^(i-1) * i * pnorm(-i*x) 582 1-(1-8*p)^k 583 }) 584 }) 585 }) 586 587 return(p) 588} 589 590 591print.efp <- function(x, ...) 592{ 593 cat("\nEmpirical Fluctuation Process:", x$type.name, "\n\n") 594 cat("Call: ") 595 print(x$call) 596 cat("\n") 597} 598 599sctest <- function(x, ...) 600{ 601 UseMethod("sctest") 602} 603 604sctest.formula <- function(formula, type = c("Rec-CUSUM", "OLS-CUSUM", 605 "Rec-MOSUM", "OLS-MOSUM", "RE", "ME", "fluctuation", "Score-CUSUM", "Nyblom-Hansen", 606 "Chow", "supF", "aveF", "expF"), h = 0.15, alt.boundary = FALSE, functional = c("max", "range", 607 "maxL2", "meanL2"), from = 0.15, to = NULL, point = 0.5, asymptotic = FALSE, data = list(), ...) 608{ 609 type <- match.arg(type) 610 if(type == "fluctuation") type <- "RE" 611 functional <- match.arg(functional) 612 dname <- paste(deparse(substitute(formula))) 613 614 if(type == "Nyblom-Hansen") { 615 type <- "Score-CUSUM" 616 functional <- "meanL2" 617 } 618 619 if(type == "Chow") 620 { 621 mf <- model.frame(formula, data = data) 622 y <- model.response(mf) 623 modelterms <- terms(formula, data = data) 624 X <- model.matrix(modelterms, data = data) 625 METHOD <- "Chow test" 626 k <- ncol(X) 627 n <- length(y) 628 629 ytsp <- NULL 630 if(is.ts(data)){ 631 if(NROW(data) == n) { 632 ytime <- time(data) 633 ytsp <- tsp(data) 634 } 635 } else { 636 env <- environment(formula) 637 if(missing(data)) data <- env 638 orig.y <- eval(attr(modelterms, "variables")[[2]], data, env) 639 if(is.ts(orig.y) & NROW(orig.y) == n){ 640 ytime <- time(orig.y) 641 ytsp <- tsp(orig.y) 642 } 643 } 644 645 ts.eps <- getOption("ts.eps") 646 647 if(length(point) > 1) { 648 if(!is.null(ytsp) && point[2] <= ytsp[3]) 649 point <- which(abs(ytime-(point[1]+(point[2]-1)/ytsp[3])) < ts.eps) 650 else 651 stop(paste(sQuote("point"), "does not specify a valid potential change point")) 652 } 653 else { 654 if(point < 1) 655 { 656 point <- floor(point*n) 657 } 658 } 659 660 if(!((point>k) & (point<(n-k)))) stop("inadmissable change point") 661 e <- lm.fit(X,y)$residuals 662 u <- c(lm.fit(as.matrix(X[(1:point),]),y[1:point])$residuals, 663 lm.fit(as.matrix(X[((point+1):n),]),y[((point+1):n)])$residuals) 664 STATISTIC <- ((sum(e^2)-sum(u^2))/k)/((sum(u^2))/(n-2*k)) 665 names(STATISTIC) <- "F" 666 if(asymptotic) 667 { 668 STATISTIC <- STATISTIC * k 669 PVAL <- 1 - pchisq(STATISTIC, k) 670 } 671 else 672 PVAL <- 1 - pf(STATISTIC, k, (n-2*k)) 673 RVAL <- list(statistic = STATISTIC, p.value = PVAL, method = METHOD, data.name = "formula") 674 class(RVAL) <- "htest" 675 } 676 else if(type %in% c("Rec-CUSUM", "OLS-CUSUM", "Rec-MOSUM", "OLS-MOSUM", "RE", "ME", "Score-CUSUM")) 677 { 678 process <- efp(formula, type = type, h = h, data = data, ...) 679 RVAL <- sctest(process, alt.boundary = alt.boundary, functional = functional) 680 } 681 else if(type %in% c("supF", "aveF", "expF")) 682 { 683 fs <- Fstats(formula, from = from, to = to, data=data, ...) 684 RVAL <- sctest(fs, type = type) 685 } 686 687 RVAL$data.name <- dname 688 return(RVAL) 689} 690 691sctest.efp <- function(x, alt.boundary = FALSE, functional = c("max", "range", "maxL2", "meanL2"), ...) 692{ 693 h <- x$par 694 type <- x$type 695 lim.process <- x$lim.process 696 functional <- match.arg(functional) 697 dname <- paste(deparse(substitute(x))) 698 METHOD <- x$type.name 699 x <- as.matrix(x$process) 700 if(lim.process %in% c("Brownian motion", "Brownian bridge")) 701 x <- x[-1, , drop = FALSE] 702 n <- nrow(x) 703 k <- ncol(x) 704 705 if(alt.boundary) { 706 if(functional == "max" & lim.process %in% c("Brownian motion", "Brownian bridge")) { 707 trim <- max(round(n * 0.001, digits = 0), 1) 708 j <- 1:n/n 709 if(lim.process == "Brownian bridge") { 710 x <- x * 1/sqrt(j * (1-j)) 711 x <- x[(trim+1):(n-trim), , drop = FALSE] 712 } else { 713 x <- x * 1/sqrt(j) 714 x <- x[trim:n, , drop = FALSE] 715 } 716 METHOD <- paste(METHOD, "with alternative boundaries") 717 } else { 718 if(lim.process == "Brownian motion") { 719 j <- 1:n/n 720 x <- x * 1/(1 + 2*j) 721 } 722 warning(paste("no alternative boundaries for", lim.process, "with", functional, "functional")) 723 } 724 } else { 725 if(lim.process == "Brownian motion") { 726 j <- 1:n/n 727 x <- x * 1/(1 + 2*j) 728 } 729 } 730 731 732 switch(functional, 733 "max" = { STAT <- max(abs(x)) }, 734 "range" = { 735 myrange <- function(y) diff(range(y)) 736 STAT <- max(apply(x, 2, myrange)) 737 METHOD <- paste(METHOD, "with range norm") 738 }, 739 "maxL2" = { 740 STAT <- max(rowSums(x^2)) 741 METHOD <- paste(METHOD, "with max L2 norm") 742 }, 743 "meanL2" = { 744 STAT <- mean(rowSums(x^2)) 745 METHOD <- paste(METHOD, "with mean L2 norm") 746 }) 747 748 switch(type, 749 "Rec-CUSUM" = { names(STAT) <- "S" }, 750 "OLS-CUSUM" = { names(STAT) <- "S0" }, 751 "Rec-MOSUM" = { names(STAT) <- "M" }, 752 "OLS-MOSUM" = { names(STAT) <- "M0" }, 753 "RE" = { names(STAT) <- "RE" }, 754 "ME" = { names(STAT) <- "ME" }, 755 names(STAT) <- "f(efp)") 756 757 PVAL <- pvalue.efp(STAT, lim.process, alt.boundary, functional = functional, h = h, k = k) 758 RVAL <- list(statistic = STAT, p.value = PVAL, method = METHOD, data.name = dname) 759 class(RVAL) <- "htest" 760 return(RVAL) 761} 762 763boundary <- function(x, ...) 764{ 765 UseMethod("boundary") 766} 767 768boundary.efp <- function(x, alpha = 0.05, alt.boundary = FALSE, functional = "max", ...) 769{ 770 h <- x$par 771 type <- x$type 772 lim.process <- x$lim.process 773 functional <- match.arg(functional, c("max", "range", "maxL2", "meanL2")) 774 proc <- as.matrix(x$process) 775 n <- nrow(proc) 776 k <- ncol(proc) 777 778 bound <- uniroot(function(y) {pvalue.efp(y, lim.process, alt.boundary = alt.boundary, 779 functional = functional, h = h, k = k) - alpha}, c(0,20))$root 780 switch(lim.process, 781 "Brownian motion" = { 782 if(functional == "max") { 783 if(alt.boundary) 784 bound <- sqrt((0:(n-1))/(n-1))*bound 785 else 786 bound <- bound + (2*bound*(0:(n-1))/(n-1)) 787 } else { 788 stop("only boundaries for Brownian motion with max functional available") 789 } 790 }, 791 "Brownian bridge" = { 792 if(!(functional == "range")) { 793 if(alt.boundary & functional == "max") 794 { 795 j <- (0:(n-1))/(n-1) 796 bound <- sqrt(j*(1-j))*bound 797 } 798 else 799 bound <- rep(bound,n) 800 } else { 801 stop("no boundaries Brownian bridge with range functional available") 802 } 803 }, 804 "Brownian motion increments" = { 805 if(functional == "max") { 806 bound <- rep(bound, n) 807 } else { 808 stop("only boundaries for Brownian motion increments with max functional available") 809 } 810 if(alt.boundary) warning("no alternative boundaries available for Brownian motion increments") 811 }, 812 "Brownian bridge increments" = { 813 if(functional == "max") { 814 bound <- rep(bound, n) 815 } else { 816 stop("only boundaries for Brownian bridge increments with max functional available") 817 } 818 if(alt.boundary) warning("no alternative boundaries available for Brownian bridge increments") 819 } 820 ) 821 822 bound <- ts(bound, end = end(x$process), frequency = frequency(x$process)) 823 return(bound) 824} 825 826lines.efp <- function(x, functional = "max", ...) 827{ 828 if(is.null(functional)) stop("lines() cannot be added for `functional = NULL'") 829 else functional <- match.arg(functional, c("max", "range", "maxL2", "meanL2")) 830 831 if(x$lim.process %in% c("Brownian bridge", "Brownian bridge increments")) { 832 z <- as.matrix(x$process) 833 k <- ncol(z) 834 835 switch(functional, 836 "max" = { 837 if(k > 1) { 838 z <- apply(abs(z), 1, max) 839 pos <- TRUE 840 } 841 }, 842 "range" = { stop("no plot available for range functional") }, 843 "maxL2" = { 844 if(x$lim.process == "Brownian bridge") { 845 z <- rowSums(z^2) 846 pos <- TRUE 847 } else { 848 stop("no test/plot available for max L2 functional") 849 } 850 }, 851 "meanL2" = { 852 if(x$lim.process == "Brownian bridge") { 853 z <- rowSums(z^2) 854 ave <- TRUE 855 pos <- TRUE 856 } else { 857 stop("no test/plot available for mean L2 functional") 858 } 859 }) 860 z <- ts(as.vector(z), start = start(x$process), frequency = frequency(x$process)) 861 } else { 862 z <- x$process 863 } 864 865 lines(z, ...) 866} 867 868