1 2R version 3.2.0 (2015-04-16) -- "Full of Ingredients" 3Copyright (C) 2015 The R Foundation for Statistical Computing 4Platform: x86_64-pc-linux-gnu (64-bit) 5 6R is free software and comes with ABSOLUTELY NO WARRANTY. 7You are welcome to redistribute it under certain conditions. 8Type 'license()' or 'licence()' for distribution details. 9 10R is a collaborative project with many contributors. 11Type 'contributors()' for more information and 12'citation()' on how to cite R or R packages in publications. 13 14Type 'demo()' for some demos, 'help()' for on-line help, or 15'help.start()' for an HTML browser interface to help. 16Type 'q()' to quit R. 17 18> pkgname <- "strucchange" 19> source(file.path(R.home("share"), "R", "examples-header.R")) 20> options(warn = 1) 21> library('strucchange') 22Loading required package: zoo 23 24Attaching package: 'zoo' 25 26The following objects are masked from 'package:base': 27 28 as.Date, as.Date.numeric 29 30Loading required package: sandwich 31> 32> base::assign(".oldSearch", base::search(), pos = 'CheckExEnv') 33> cleanEx() 34> nameEx("BostonHomicide") 35> ### * BostonHomicide 36> 37> flush(stderr()); flush(stdout()) 38> 39> ### Name: BostonHomicide 40> ### Title: Youth Homicides in Boston 41> ### Aliases: BostonHomicide 42> ### Keywords: datasets 43> 44> ### ** Examples 45> 46> data("BostonHomicide") 47> attach(BostonHomicide) 48> 49> ## data from Table 1 50> tapply(homicides, year, mean) 51 1992 1993 1994 1995 1996 1997 1998 523.083333 4.000000 3.166667 3.833333 2.083333 1.250000 0.800000 53> populationBM[0:6*12 + 7] 54[1] 12977 12455 12272 12222 11895 12038 NA 55> tapply(ahomicides25, year, mean) 56 1992 1993 1994 1995 1996 1997 1998 573.250000 4.166667 3.916667 4.166667 2.666667 2.333333 1.400000 58> tapply(ahomicides35, year, mean) 59 1992 1993 1994 1995 1996 1997 1998 600.8333333 1.0833333 1.3333333 1.1666667 1.0833333 0.7500000 0.4000000 61> population[0:6*12 + 7] 62[1] 228465 227218 226611 231367 230744 228696 NA 63> unemploy[0:6*12 + 7] 64[1] 20.2 18.8 15.9 14.7 13.8 12.6 NA 65> 66> ## model A 67> ## via OLS 68> fmA <- lm(homicides ~ populationBM + season) 69> anova(fmA) 70Analysis of Variance Table 71 72Response: homicides 73 Df Sum Sq Mean Sq F value Pr(>F) 74populationBM 1 14.364 14.3642 3.7961 0.05576 . 75season 11 47.254 4.2959 1.1353 0.34985 76Residuals 64 242.174 3.7840 77--- 78Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 79> ## as GLM 80> fmA1 <- glm(homicides ~ populationBM + season, family = poisson) 81> anova(fmA1, test = "Chisq") 82Analysis of Deviance Table 83 84Model: poisson, link: log 85 86Response: homicides 87 88Terms added sequentially (first to last) 89 90 91 Df Deviance Resid. Df Resid. Dev Pr(>Chi) 92NULL 76 115.649 93populationBM 1 4.9916 75 110.657 0.02547 * 94season 11 18.2135 64 92.444 0.07676 . 95--- 96Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 97> 98> ## model B & C 99> fmB <- lm(homicides ~ populationBM + season + ahomicides25) 100> fmC <- lm(homicides ~ populationBM + season + ahomicides25 + unemploy) 101> 102> detach(BostonHomicide) 103> 104> 105> 106> cleanEx() 107> nameEx("DJIA") 108> ### * DJIA 109> 110> flush(stderr()); flush(stdout()) 111> 112> ### Name: DJIA 113> ### Title: Dow Jones Industrial Average 114> ### Aliases: DJIA 115> ### Keywords: datasets 116> 117> ### ** Examples 118> 119> data("DJIA") 120> ## look at log-difference returns 121> djia <- diff(log(DJIA)) 122> plot(djia) 123> 124> ## convenience functions 125> ## set up a normal regression model which 126> ## explicitely also models the variance 127> normlm <- function(formula, data = list()) { 128+ rval <- lm(formula, data = data) 129+ class(rval) <- c("normlm", "lm") 130+ return(rval) 131+ } 132> estfun.normlm <- function(obj) { 133+ res <- residuals(obj) 134+ ef <- NextMethod(obj) 135+ sigma2 <- mean(res^2) 136+ rval <- cbind(ef, res^2 - sigma2) 137+ colnames(rval) <- c(colnames(ef), "(Variance)") 138+ return(rval) 139+ } 140> 141> ## normal model (with constant mean and variance) for log returns 142> m1 <- gefp(djia ~ 1, fit = normlm, vcov = meatHAC, sandwich = FALSE) 143> plot(m1, aggregate = FALSE) 144> ## suggests a clear break in the variance (but not the mean) 145> 146> ## dating 147> bp <- breakpoints(I(djia^2) ~ 1) 148> plot(bp) 149> ## -> clearly one break 150> bp 151 152 Optimal 2-segment partition: 153 154Call: 155breakpoints.formula(formula = I(djia^2) ~ 1) 156 157Breakpoints at observation number: 15889 159 160Corresponding to breakdates: 1610.552795 162> time(djia)[bp$breakpoints] 163[1] "1973-03-16" 164> 165> ## visualization 166> plot(djia) 167> abline(v = time(djia)[bp$breakpoints], lty = 2) 168> lines(time(djia)[confint(bp)$confint[c(1,3)]], rep(min(djia), 2), col = 2, type = "b", pch = 3) 169> 170> 171> 172> cleanEx() 173> nameEx("Fstats") 174> ### * Fstats 175> 176> flush(stderr()); flush(stdout()) 177> 178> ### Name: Fstats 179> ### Title: F Statistics 180> ### Aliases: Fstats print.Fstats 181> ### Keywords: regression 182> 183> ### ** Examples 184> 185> ## Nile data with one breakpoint: the annual flows drop in 1898 186> ## because the first Ashwan dam was built 187> data("Nile") 188> plot(Nile) 189> 190> ## test the null hypothesis that the annual flow remains constant 191> ## over the years 192> fs.nile <- Fstats(Nile ~ 1) 193> plot(fs.nile) 194> sctest(fs.nile) 195 196 supF test 197 198data: fs.nile 199sup.F = 75.93, p-value = 2.22e-16 200 201> ## visualize the breakpoint implied by the argmax of the F statistics 202> plot(Nile) 203> lines(breakpoints(fs.nile)) 204> 205> ## UK Seatbelt data: a SARIMA(1,0,0)(1,0,0)_12 model 206> ## (fitted by OLS) is used and reveals (at least) two 207> ## breakpoints - one in 1973 associated with the oil crisis and 208> ## one in 1983 due to the introduction of compulsory 209> ## wearing of seatbelts in the UK. 210> data("UKDriverDeaths") 211> seatbelt <- log10(UKDriverDeaths) 212> seatbelt <- cbind(seatbelt, lag(seatbelt, k = -1), lag(seatbelt, k = -12)) 213> colnames(seatbelt) <- c("y", "ylag1", "ylag12") 214> seatbelt <- window(seatbelt, start = c(1970, 1), end = c(1984,12)) 215> plot(seatbelt[,"y"], ylab = expression(log[10](casualties))) 216> 217> ## compute F statistics for potential breakpoints between 218> ## 1971(6) (corresponds to from = 0.1) and 1983(6) (corresponds to 219> ## to = 0.9 = 1 - from, the default) 220> ## compute F statistics 221> fs <- Fstats(y ~ ylag1 + ylag12, data = seatbelt, from = 0.1) 222> ## this gives the same result 223> fs <- Fstats(y ~ ylag1 + ylag12, data = seatbelt, from = c(1971, 6), 224+ to = c(1983, 6)) 225> ## plot the F statistics 226> plot(fs, alpha = 0.01) 227> ## plot F statistics with aveF boundary 228> plot(fs, aveF = TRUE) 229> ## perform the expF test 230> sctest(fs, type = "expF") 231 232 expF test 233 234data: fs 235exp.F = 6.4247, p-value = 0.008093 236 237> 238> 239> 240> cleanEx() 241> nameEx("GermanM1") 242> ### * GermanM1 243> 244> flush(stderr()); flush(stdout()) 245> 246> ### Encoding: UTF-8 247> 248> ### Name: GermanM1 249> ### Title: German M1 Money Demand 250> ### Aliases: GermanM1 historyM1 monitorM1 251> ### Keywords: datasets 252> 253> ### ** Examples 254> 255> data("GermanM1") 256> ## Lütkepohl et al. (1999) use the following model 257> LTW.model <- dm ~ dy2 + dR + dR1 + dp + m1 + y1 + R1 + season 258> ## Zeileis et al. (2005) use 259> M1.model <- dm ~ dy2 + dR + dR1 + dp + ecm.res + season 260> 261> 262> ## historical tests 263> ols <- efp(LTW.model, data = GermanM1, type = "OLS-CUSUM") 264> plot(ols) 265> re <- efp(LTW.model, data = GermanM1, type = "fluctuation") 266> plot(re) 267> fs <- Fstats(LTW.model, data = GermanM1, from = 0.1) 268> plot(fs) 269> 270> ## monitoring 271> M1 <- historyM1 272> ols.efp <- efp(M1.model, type = "OLS-CUSUM", data = M1) 273> newborder <- function(k) 1.5778*k/118 274> ols.mefp <- mefp(ols.efp, period = 2) 275> ols.mefp2 <- mefp(ols.efp, border = newborder) 276> M1 <- GermanM1 277> ols.mon <- monitor(ols.mefp) 278Break detected at observation # 128 279> ols.mon2 <- monitor(ols.mefp2) 280Break detected at observation # 135 281> plot(ols.mon) 282> lines(boundary(ols.mon2), col = 2) 283> 284> ## dating 285> bp <- breakpoints(LTW.model, data = GermanM1) 286> summary(bp) 287 288 Optimal (m+1)-segment partition: 289 290Call: 291breakpoints.formula(formula = LTW.model, data = GermanM1) 292 293Breakpoints at observation number: 294 295m = 1 119 296m = 2 42 119 297m = 3 48 71 119 298m = 4 27 48 71 119 299m = 5 27 48 71 98 119 300 301Corresponding to breakdates: 302 303m = 1 1990(3) 304m = 2 1971(2) 1990(3) 305m = 3 1972(4) 1978(3) 1990(3) 306m = 4 1967(3) 1972(4) 1978(3) 1990(3) 307m = 5 1967(3) 1972(4) 1978(3) 1985(2) 1990(3) 308 309Fit: 310 311m 0 1 2 3 4 5 312RSS 3.683e-02 1.916e-02 1.522e-02 1.301e-02 1.053e-02 9.198e-03 313BIC -6.974e+02 -7.296e+02 -7.025e+02 -6.653e+02 -6.356e+02 -5.952e+02 314> plot(bp) 315> 316> plot(fs) 317> lines(confint(bp)) 318> 319> 320> 321> cleanEx() 322> nameEx("Grossarl") 323> ### * Grossarl 324> 325> flush(stderr()); flush(stdout()) 326> 327> ### Name: Grossarl 328> ### Title: Marriages, Births and Deaths in Grossarl 329> ### Aliases: Grossarl 330> ### Keywords: datasets 331> 332> ### ** Examples 333> 334> data("Grossarl") 335> 336> ## time series of births, deaths, marriages 337> ########################################### 338> 339> with(Grossarl, plot(cbind(deaths, illegitimate + legitimate, marriages), 340+ plot.type = "single", col = grey(c(0.7, 0, 0)), lty = c(1, 1, 3), 341+ lwd = 1.5, ylab = "annual Grossarl series")) 342> legend("topright", c("deaths", "births", "marriages"), col = grey(c(0.7, 0, 0)), 343+ lty = c(1, 1, 3), bty = "n") 344> 345> ## illegitimate births 346> ###################### 347> ## lm + MOSUM 348> plot(Grossarl$fraction) 349> fm.min <- lm(fraction ~ politics, data = Grossarl) 350> fm.ext <- lm(fraction ~ politics + morals + nuptiality + marriages, 351+ data = Grossarl) 352> lines(ts(fitted(fm.min), start = 1700), col = 2) 353> lines(ts(fitted(fm.ext), start = 1700), col = 4) 354> mos.min <- efp(fraction ~ politics, data = Grossarl, type = "OLS-MOSUM") 355> mos.ext <- efp(fraction ~ politics + morals + nuptiality + marriages, 356+ data = Grossarl, type = "OLS-MOSUM") 357> plot(mos.min) 358> lines(mos.ext, lty = 2) 359> 360> ## dating 361> bp <- breakpoints(fraction ~ 1, data = Grossarl, h = 0.1) 362> summary(bp) 363 364 Optimal (m+1)-segment partition: 365 366Call: 367breakpoints.formula(formula = fraction ~ 1, h = 0.1, data = Grossarl) 368 369Breakpoints at observation number: 370 371m = 1 127 372m = 2 55 122 373m = 3 55 124 180 374m = 4 55 122 157 179 375m = 5 54 86 122 157 179 376m = 6 35 55 86 122 157 179 377m = 7 35 55 80 101 122 157 179 378m = 8 35 55 79 99 119 139 159 179 379 380Corresponding to breakdates: 381 382m = 1 1826 383m = 2 1754 1821 384m = 3 1754 1823 1879 385m = 4 1754 1821 1856 1878 386m = 5 1753 1785 1821 1856 1878 387m = 6 1734 1754 1785 1821 1856 1878 388m = 7 1734 1754 1779 1800 1821 1856 1878 389m = 8 1734 1754 1778 1798 1818 1838 1858 1878 390 391Fit: 392 393m 0 1 2 3 4 5 6 394RSS 1.1088 0.8756 0.6854 0.6587 0.6279 0.6019 0.5917 395BIC -460.8402 -497.4625 -535.8459 -533.1857 -532.1789 -530.0501 -522.8510 396 397m 7 8 398RSS 0.5934 0.6084 399BIC -511.7017 -496.0924 400> ## RSS, BIC, AIC 401> plot(bp) 402> plot(0:8, AIC(bp), type = "b") 403> 404> ## probably use 5 or 6 breakpoints and compare with 405> ## coding of the factors as used by us 406> ## 407> ## politics 1803 1816 1850 408> ## morals 1736 1753 1771 1803 409> ## nuptiality 1803 1810 1816 1883 410> ## 411> ## m = 5 1753 1785 1821 1856 1878 412> ## m = 6 1734 1754 1785 1821 1856 1878 413> ## 6 2 5 1 4 3 414> 415> ## fitted models 416> coef(bp, breaks = 6) 417 (Intercept) 4181700 - 1734 0.16933985 4191735 - 1754 0.14078070 4201755 - 1785 0.09890276 4211786 - 1821 0.05955620 4221822 - 1856 0.17441529 4231857 - 1878 0.22425604 4241879 - 1899 0.15414723 425> plot(Grossarl$fraction) 426> lines(fitted(bp, breaks = 6), col = 2) 427> lines(ts(fitted(fm.ext), start = 1700), col = 4) 428> 429> 430> ## marriages 431> ############ 432> ## lm + MOSUM 433> plot(Grossarl$marriages) 434> fm.min <- lm(marriages ~ politics, data = Grossarl) 435> fm.ext <- lm(marriages ~ politics + morals + nuptiality, data = Grossarl) 436> lines(ts(fitted(fm.min), start = 1700), col = 2) 437> lines(ts(fitted(fm.ext), start = 1700), col = 4) 438> mos.min <- efp(marriages ~ politics, data = Grossarl, type = "OLS-MOSUM") 439> mos.ext <- efp(marriages ~ politics + morals + nuptiality, data = Grossarl, 440+ type = "OLS-MOSUM") 441> plot(mos.min) 442> lines(mos.ext, lty = 2) 443> 444> ## dating 445> bp <- breakpoints(marriages ~ 1, data = Grossarl, h = 0.1) 446> summary(bp) 447 448 Optimal (m+1)-segment partition: 449 450Call: 451breakpoints.formula(formula = marriages ~ 1, h = 0.1, data = Grossarl) 452 453Breakpoints at observation number: 454 455m = 1 114 456m = 2 39 114 457m = 3 39 114 176 458m = 4 39 95 115 176 459m = 5 39 62 95 115 176 460m = 6 39 62 95 115 136 176 461m = 7 39 62 95 115 136 156 176 462m = 8 21 41 62 95 115 136 156 176 463 464Corresponding to breakdates: 465 466m = 1 1813 467m = 2 1738 1813 468m = 3 1738 1813 1875 469m = 4 1738 1794 1814 1875 470m = 5 1738 1761 1794 1814 1875 471m = 6 1738 1761 1794 1814 1835 1875 472m = 7 1738 1761 1794 1814 1835 1855 1875 473m = 8 1720 1740 1761 1794 1814 1835 1855 1875 474 475Fit: 476 477m 0 1 2 3 4 5 6 7 8 478RSS 3832 3059 2863 2723 2671 2634 2626 2626 2645 479BIC 1169 1134 1132 1132 1139 1147 1157 1167 1179 480> ## RSS, BIC, AIC 481> plot(bp) 482> plot(0:8, AIC(bp), type = "b") 483> 484> ## probably use 3 or 4 breakpoints and compare with 485> ## coding of the factors as used by us 486> ## 487> ## politics 1803 1816 1850 488> ## morals 1736 1753 1771 1803 489> ## nuptiality 1803 1810 1816 1883 490> ## 491> ## m = 3 1738 1813 1875 492> ## m = 4 1738 1794 1814 1875 493> ## 2 4 1 3 494> 495> ## fitted models 496> coef(bp, breaks = 4) 497 (Intercept) 4981700 - 1738 13.487179 4991739 - 1794 10.160714 5001795 - 1814 12.150000 5011815 - 1875 6.885246 5021876 - 1899 9.750000 503> plot(Grossarl$marriages) 504> lines(fitted(bp, breaks = 4), col = 2) 505> lines(ts(fitted(fm.ext), start = 1700), col = 4) 506> 507> 508> 509> cleanEx() 510> nameEx("PhillipsCurve") 511> ### * PhillipsCurve 512> 513> flush(stderr()); flush(stdout()) 514> 515> ### Name: PhillipsCurve 516> ### Title: UK Phillips Curve Equation Data 517> ### Aliases: PhillipsCurve 518> ### Keywords: datasets 519> 520> ### ** Examples 521> 522> ## load and plot data 523> data("PhillipsCurve") 524> uk <- window(PhillipsCurve, start = 1948) 525> plot(uk[, "dp"]) 526> 527> ## AR(1) inflation model 528> ## estimate breakpoints 529> bp.inf <- breakpoints(dp ~ dp1, data = uk, h = 8) 530> plot(bp.inf) 531> summary(bp.inf) 532 533 Optimal (m+1)-segment partition: 534 535Call: 536breakpoints.formula(formula = dp ~ dp1, h = 8, data = uk) 537 538Breakpoints at observation number: 539 540m = 1 20 541m = 2 20 28 542m = 3 9 20 28 543 544Corresponding to breakdates: 545 546m = 1 1967 547m = 2 1967 1975 548m = 3 1956 1967 1975 549 550Fit: 551 552m 0 1 2 3 553RSS 0.03068 0.02672 0.01838 0.01786 554BIC -162.34174 -156.80265 -160.70385 -150.78479 555> 556> ## fit segmented model with three breaks 557> fac.inf <- breakfactor(bp.inf, breaks = 2, label = "seg") 558> fm.inf <- lm(dp ~ 0 + fac.inf/dp1, data = uk) 559> summary(fm.inf) 560 561Call: 562lm(formula = dp ~ 0 + fac.inf/dp1, data = uk) 563 564Residuals: 565 Min 1Q Median 3Q Max 566-0.046987 -0.014861 -0.003593 0.006286 0.058081 567 568Coefficients: 569 Estimate Std. Error t value Pr(>|t|) 570fac.infseg1 0.024501 0.011176 2.192 0.0353 * 571fac.infseg2 -0.000775 0.017853 -0.043 0.9656 572fac.infseg3 0.017603 0.015007 1.173 0.2489 573fac.infseg1:dp1 0.274012 0.269892 1.015 0.3171 574fac.infseg2:dp1 1.343369 0.224521 5.983 9.05e-07 *** 575fac.infseg3:dp1 0.683410 0.130106 5.253 8.07e-06 *** 576--- 577Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 578 579Residual standard error: 0.02325 on 34 degrees of freedom 580Multiple R-squared: 0.9237, Adjusted R-squared: 0.9103 581F-statistic: 68.64 on 6 and 34 DF, p-value: < 2.2e-16 582 583> 584> ## Results from Table 2 in Bai & Perron (2003): 585> ## coefficient estimates 586> coef(bp.inf, breaks = 2) 587 (Intercept) dp1 5881948 - 1967 0.0245010729 0.2740125 5891968 - 1975 -0.0007750299 1.3433686 5901976 - 1987 0.0176032179 0.6834098 591> ## corresponding standard errors 592> sqrt(sapply(vcov(bp.inf, breaks = 2), diag)) 593 1948 - 1967 1968 - 1975 1976 - 1987 594(Intercept) 0.008268814 0.01985539 0.01571339 595dp1 0.199691273 0.24969992 0.13622996 596> ## breakpoints and confidence intervals 597> confint(bp.inf, breaks = 2) 598 599 Confidence intervals for breakpoints 600 of optimal 3-segment partition: 601 602Call: 603confint.breakpointsfull(object = bp.inf, breaks = 2) 604 605Breakpoints at observation number: 606 2.5 % breakpoints 97.5 % 6071 18 20 25 6082 26 28 34 609 610Corresponding to breakdates: 611 2.5 % breakpoints 97.5 % 6121 1965 1967 1972 6132 1973 1975 1981 614> 615> ## Phillips curve equation 616> ## estimate breakpoints 617> bp.pc <- breakpoints(dw ~ dp1 + du + u1, data = uk, h = 5, breaks = 5) 618> ## look at RSS and BIC 619> plot(bp.pc) 620> summary(bp.pc) 621 622 Optimal (m+1)-segment partition: 623 624Call: 625breakpoints.formula(formula = dw ~ dp1 + du + u1, h = 5, breaks = 5, 626 data = uk) 627 628Breakpoints at observation number: 629 630m = 1 26 631m = 2 20 28 632m = 3 9 25 30 633m = 4 11 16 25 30 634m = 5 11 16 22 27 32 635 636Corresponding to breakdates: 637 638m = 1 1973 639m = 2 1967 1975 640m = 3 1956 1972 1977 641m = 4 1958 1963 1972 1977 642m = 5 1958 1963 1969 1974 1979 643 644Fit: 645 646m 0 1 2 3 4 5 647RSS 3.409e-02 1.690e-02 1.062e-02 7.835e-03 5.183e-03 3.388e-03 648BIC -1.508e+02 -1.604e+02 -1.605e+02 -1.542e+02 -1.523e+02 -1.509e+02 649> 650> ## fit segmented model with three breaks 651> fac.pc <- breakfactor(bp.pc, breaks = 2, label = "seg") 652> fm.pc <- lm(dw ~ 0 + fac.pc/dp1 + du + u1, data = uk) 653> summary(fm.pc) 654 655Call: 656lm(formula = dw ~ 0 + fac.pc/dp1 + du + u1, data = uk) 657 658Residuals: 659 Min 1Q Median 3Q Max 660-0.041392 -0.011516 0.000089 0.010036 0.044539 661 662Coefficients: 663 Estimate Std. Error t value Pr(>|t|) 664fac.pcseg1 0.06574 0.01169 5.623 3.24e-06 *** 665fac.pcseg2 0.06231 0.01883 3.310 0.00232 ** 666fac.pcseg3 0.18093 0.05388 3.358 0.00204 ** 667du -0.14408 0.58218 -0.247 0.80611 668u1 -0.87516 0.37274 -2.348 0.02523 * 669fac.pcseg1:dp1 0.09373 0.24053 0.390 0.69936 670fac.pcseg2:dp1 1.23143 0.20498 6.008 1.06e-06 *** 671fac.pcseg3:dp1 0.01618 0.25667 0.063 0.95013 672--- 673Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 674 675Residual standard error: 0.02021 on 32 degrees of freedom 676Multiple R-squared: 0.9655, Adjusted R-squared: 0.9569 677F-statistic: 112 on 8 and 32 DF, p-value: < 2.2e-16 678 679> 680> ## Results from Table 3 in Bai & Perron (2003): 681> ## coefficient estimates 682> coef(fm.pc) 683 fac.pcseg1 fac.pcseg2 fac.pcseg3 du u1 684 0.06574278 0.06231337 0.18092502 -0.14408073 -0.87515585 685fac.pcseg1:dp1 fac.pcseg2:dp1 fac.pcseg3:dp1 686 0.09372759 1.23143008 0.01617826 687> ## corresponding standard errors 688> sqrt(diag(vcov(fm.pc))) 689 fac.pcseg1 fac.pcseg2 fac.pcseg3 du u1 690 0.01169149 0.01882668 0.05388166 0.58217571 0.37273955 691fac.pcseg1:dp1 fac.pcseg2:dp1 fac.pcseg3:dp1 692 0.24052539 0.20497973 0.25666903 693> ## breakpoints and confidence intervals 694> confint(bp.pc, breaks = 2, het.err = FALSE) 695 696 Confidence intervals for breakpoints 697 of optimal 3-segment partition: 698 699Call: 700confint.breakpointsfull(object = bp.pc, breaks = 2, het.err = FALSE) 701 702Breakpoints at observation number: 703 2.5 % breakpoints 97.5 % 7041 19 20 21 7052 27 28 29 706 707Corresponding to breakdates: 708 2.5 % breakpoints 97.5 % 7091 1966 1967 1968 7102 1974 1975 1976 711> 712> 713> 714> cleanEx() 715> nameEx("RealInt") 716> ### * RealInt 717> 718> flush(stderr()); flush(stdout()) 719> 720> ### Name: RealInt 721> ### Title: US Ex-post Real Interest Rate 722> ### Aliases: RealInt 723> ### Keywords: datasets 724> 725> ### ** Examples 726> 727> ## load and plot data 728> data("RealInt") 729> plot(RealInt) 730> 731> ## estimate breakpoints 732> bp.ri <- breakpoints(RealInt ~ 1, h = 15) 733> plot(bp.ri) 734> summary(bp.ri) 735 736 Optimal (m+1)-segment partition: 737 738Call: 739breakpoints.formula(formula = RealInt ~ 1, h = 15) 740 741Breakpoints at observation number: 742 743m = 1 79 744m = 2 47 79 745m = 3 24 47 79 746m = 4 24 47 64 79 747m = 5 16 31 47 64 79 748 749Corresponding to breakdates: 750 751m = 1 1980(3) 752m = 2 1972(3) 1980(3) 753m = 3 1966(4) 1972(3) 1980(3) 754m = 4 1966(4) 1972(3) 1976(4) 1980(3) 755m = 5 1964(4) 1968(3) 1972(3) 1976(4) 1980(3) 756 757Fit: 758 759m 0 1 2 3 4 5 760RSS 1214.9 645.0 456.0 445.2 444.9 449.6 761BIC 555.7 499.8 473.3 480.1 489.3 499.7 762> 763> ## fit segmented model with three breaks 764> fac.ri <- breakfactor(bp.ri, breaks = 3, label = "seg") 765> fm.ri <- lm(RealInt ~ 0 + fac.ri) 766> summary(fm.ri) 767 768Call: 769lm(formula = RealInt ~ 0 + fac.ri) 770 771Residuals: 772 Min 1Q Median 3Q Max 773-4.5157 -1.3674 -0.0578 1.3248 6.0990 774 775Coefficients: 776 Estimate Std. Error t value Pr(>|t|) 777fac.riseg1 1.8236 0.4329 4.213 5.57e-05 *** 778fac.riseg2 0.8661 0.4422 1.959 0.053 . 779fac.riseg3 -1.7961 0.3749 -4.791 5.83e-06 *** 780fac.riseg4 5.6429 0.4329 13.036 < 2e-16 *** 781--- 782Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 783 784Residual standard error: 2.121 on 99 degrees of freedom 785Multiple R-squared: 0.6842, Adjusted R-squared: 0.6714 786F-statistic: 53.62 on 4 and 99 DF, p-value: < 2.2e-16 787 788> 789> ## setup kernel HAC estimator 790> vcov.ri <- function(x, ...) kernHAC(x, kernel = "Quadratic Spectral", 791+ prewhite = 1, approx = "AR(1)", ...) 792> 793> ## Results from Table 1 in Bai & Perron (2003): 794> ## coefficient estimates 795> coef(bp.ri, breaks = 3) 796 (Intercept) 7971961(1) - 1966(4) 1.8236167 7981967(1) - 1972(3) 0.8660848 7991972(4) - 1980(3) -1.7961384 8001980(4) - 1986(3) 5.6428896 801> ## corresponding standard errors 802> sapply(vcov(bp.ri, breaks = 3, vcov = vcov.ri), sqrt) 8031961(1) - 1966(4) 1967(1) - 1972(3) 1972(4) - 1980(3) 1980(4) - 1986(3) 804 0.1857577 0.1499849 0.5026749 0.5887460 805> ## breakpoints and confidence intervals 806> confint(bp.ri, breaks = 3, vcov = vcov.ri) 807 808 Confidence intervals for breakpoints 809 of optimal 4-segment partition: 810 811Call: 812confint.breakpointsfull(object = bp.ri, breaks = 3, vcov. = vcov.ri) 813 814Breakpoints at observation number: 815 2.5 % breakpoints 97.5 % 8161 18 24 35 8172 33 47 48 8183 77 79 81 819 820Corresponding to breakdates: 821Warning: Overlapping confidence intervals 822 2.5 % breakpoints 97.5 % 8231 1965(2) 1966(4) 1969(3) 8242 1969(1) 1972(3) 1972(4) 8253 1980(1) 1980(3) 1981(1) 826> 827> ## Visualization 828> plot(RealInt) 829> lines(as.vector(time(RealInt)), fitted(fm.ri), col = 4) 830> lines(confint(bp.ri, breaks = 3, vcov = vcov.ri)) 831Warning: Overlapping confidence intervals 832> 833> 834> 835> cleanEx() 836> nameEx("SP2001") 837> ### * SP2001 838> 839> flush(stderr()); flush(stdout()) 840> 841> ### Name: SP2001 842> ### Title: S&P 500 Stock Prices 843> ### Aliases: SP2001 844> ### Keywords: datasets 845> 846> ### ** Examples 847> 848> ## load and transform data 849> ## (DAL: Delta Air Lines, LU: Lucent Technologies) 850> data("SP2001") 851> stock.prices <- SP2001[, c("DAL", "LU")] 852> stock.returns <- diff(log(stock.prices)) 853> 854> ## price and return series 855> plot(stock.prices, ylab = c("Delta Air Lines", "Lucent Technologies"), main = "") 856> plot(stock.returns, ylab = c("Delta Air Lines", "Lucent Technologies"), main = "") 857> 858> ## monitoring of DAL series 859> myborder <- function(k) 1.939*k/28 860> x <- as.vector(stock.returns[, "DAL"][1:28]) 861> dal.cusum <- mefp(x ~ 1, type = "OLS-CUSUM", border = myborder) 862> dal.mosum <- mefp(x ~ 1, type = "OLS-MOSUM", h = 0.5, period = 4) 863> x <- as.vector(stock.returns[, "DAL"]) 864> dal.cusum <- monitor(dal.cusum) 865Break detected at observation # 29 866> dal.mosum <- monitor(dal.mosum) 867Break detected at observation # 29 868> 869> ## monitoring of LU series 870> x <- as.vector(stock.returns[, "LU"][1:28]) 871> lu.cusum <- mefp(x ~ 1, type = "OLS-CUSUM", border = myborder) 872> lu.mosum <- mefp(x ~ 1, type = "OLS-MOSUM", h = 0.5, period = 4) 873> x <- as.vector(stock.returns[, "LU"]) 874> lu.cusum <- monitor(lu.cusum) 875> lu.mosum <- monitor(lu.mosum) 876> 877> ## pretty plotting 878> ## (needs some work because lm() does not keep "zoo" attributes) 879> cus.bound <- zoo(c(rep(NA, 27), myborder(28:102)), index(stock.returns)) 880> mos.bound <- as.vector(boundary(dal.mosum)) 881> mos.bound <- zoo(c(rep(NA, 27), mos.bound[1], mos.bound), index(stock.returns)) 882> 883> ## Lucent Technologies: CUSUM test 884> plot(zoo(c(lu.cusum$efpprocess, lu.cusum$process), index(stock.prices)), 885+ ylim = c(-1, 1) * coredata(cus.bound)[102], xlab = "Time", ylab = "empirical fluctuation process") 886> abline(0, 0) 887> abline(v = as.Date("2001-09-10"), lty = 2) 888> lines(cus.bound, col = 2) 889> lines(-cus.bound, col = 2) 890> 891> ## Lucent Technologies: MOSUM test 892> plot(zoo(c(lu.mosum$efpprocess, lu.mosum$process), index(stock.prices)[-(1:14)]), 893+ ylim = c(-1, 1) * coredata(mos.bound)[102], xlab = "Time", ylab = "empirical fluctuation process") 894> abline(0, 0) 895> abline(v = as.Date("2001-09-10"), lty = 2) 896> lines(mos.bound, col = 2) 897> lines(-mos.bound, col = 2) 898> 899> ## Delta Air Lines: CUSUM test 900> plot(zoo(c(dal.cusum$efpprocess, dal.cusum$process), index(stock.prices)), 901+ ylim = c(-1, 1) * coredata(cus.bound)[102], xlab = "Time", ylab = "empirical fluctuation process") 902> abline(0, 0) 903> abline(v = as.Date("2001-09-10"), lty = 2) 904> lines(cus.bound, col = 2) 905> lines(-cus.bound, col = 2) 906> 907> ## Delta Air Lines: MOSUM test 908> plot(zoo(c(dal.mosum$efpprocess, dal.mosum$process), index(stock.prices)[-(1:14)]), 909+ ylim = range(dal.mosum$process), xlab = "Time", ylab = "empirical fluctuation process") 910> abline(0, 0) 911> abline(v = as.Date("2001-09-10"), lty = 2) 912> lines(mos.bound, col = 2) 913> lines(-mos.bound, col = 2) 914> 915> 916> 917> cleanEx() 918> nameEx("USIncExp") 919> ### * USIncExp 920> 921> flush(stderr()); flush(stdout()) 922> 923> ### Name: USIncExp 924> ### Title: Income and Expenditures in the US 925> ### Aliases: USIncExp 926> ### Keywords: datasets 927> 928> ### ** Examples 929> 930> ## These example are presented in the vignette distributed with this 931> ## package, the code was generated by Stangle("strucchange-intro.Rnw") 932> 933> ################################################### 934> ### chunk number 1: data 935> ################################################### 936> library("strucchange") 937> data("USIncExp") 938> plot(USIncExp, plot.type = "single", col = 1:2, ylab = "billion US$") 939> legend(1960, max(USIncExp), c("income", "expenditures"), 940+ lty = c(1,1), col = 1:2, bty = "n") 941> 942> 943> ################################################### 944> ### chunk number 2: subset 945> ################################################### 946> library("strucchange") 947> data("USIncExp") 948> USIncExp2 <- window(USIncExp, start = c(1985,12)) 949> 950> 951> ################################################### 952> ### chunk number 3: ecm-setup 953> ################################################### 954> coint.res <- residuals(lm(expenditure ~ income, data = USIncExp2)) 955> coint.res <- lag(ts(coint.res, start = c(1985,12), freq = 12), k = -1) 956> USIncExp2 <- cbind(USIncExp2, diff(USIncExp2), coint.res) 957> USIncExp2 <- window(USIncExp2, start = c(1986,1), end = c(2001,2)) 958> colnames(USIncExp2) <- c("income", "expenditure", "diff.income", 959+ "diff.expenditure", "coint.res") 960> ecm.model <- diff.expenditure ~ coint.res + diff.income 961> 962> 963> ################################################### 964> ### chunk number 4: ts-used 965> ################################################### 966> plot(USIncExp2[,3:5], main = "") 967> 968> 969> ################################################### 970> ### chunk number 5: efp 971> ################################################### 972> ocus <- efp(ecm.model, type="OLS-CUSUM", data=USIncExp2) 973> me <- efp(ecm.model, type="ME", data=USIncExp2, h=0.2) 974> 975> 976> ################################################### 977> ### chunk number 6: efp-boundary 978> ################################################### 979> bound.ocus <- boundary(ocus, alpha=0.05) 980> 981> 982> ################################################### 983> ### chunk number 7: OLS-CUSUM 984> ################################################### 985> plot(ocus) 986> 987> 988> ################################################### 989> ### chunk number 8: efp-boundary2 990> ################################################### 991> plot(ocus, boundary = FALSE) 992> lines(bound.ocus, col = 4) 993> lines(-bound.ocus, col = 4) 994> 995> 996> ################################################### 997> ### chunk number 9: ME-null 998> ################################################### 999> plot(me, functional = NULL) 1000> 1001> 1002> ################################################### 1003> ### chunk number 10: efp-sctest 1004> ################################################### 1005> sctest(ocus) 1006 1007 OLS-based CUSUM test 1008 1009data: ocus 1010S0 = 1.5511, p-value = 0.01626 1011 1012> 1013> 1014> ################################################### 1015> ### chunk number 11: efp-sctest2 1016> ################################################### 1017> sctest(ecm.model, type="OLS-CUSUM", data=USIncExp2) 1018 1019 OLS-based CUSUM test 1020 1021data: ecm.model 1022S0 = 1.5511, p-value = 0.01626 1023 1024> 1025> 1026> ################################################### 1027> ### chunk number 12: Fstats 1028> ################################################### 1029> fs <- Fstats(ecm.model, from = c(1990, 1), to = c(1999,6), data = USIncExp2) 1030> 1031> 1032> ################################################### 1033> ### chunk number 13: Fstats-plot 1034> ################################################### 1035> plot(fs) 1036> 1037> 1038> ################################################### 1039> ### chunk number 14: pval-plot 1040> ################################################### 1041> plot(fs, pval=TRUE) 1042> 1043> 1044> ################################################### 1045> ### chunk number 15: aveF-plot 1046> ################################################### 1047> plot(fs, aveF=TRUE) 1048> 1049> 1050> ################################################### 1051> ### chunk number 16: Fstats-sctest 1052> ################################################### 1053> sctest(fs, type="expF") 1054 1055 expF test 1056 1057data: fs 1058exp.F = 8.9955, p-value = 0.001311 1059 1060> 1061> 1062> ################################################### 1063> ### chunk number 17: Fstats-sctest2 1064> ################################################### 1065> sctest(ecm.model, type = "expF", from = 49, to = 162, data = USIncExp2) 1066 1067 expF test 1068 1069data: ecm.model 1070exp.F = 8.9955, p-value = 0.001311 1071 1072> 1073> 1074> ################################################### 1075> ### chunk number 18: mefp 1076> ################################################### 1077> USIncExp3 <- window(USIncExp2, start = c(1986, 1), end = c(1989,12)) 1078> me.mefp <- mefp(ecm.model, type = "ME", data = USIncExp3, alpha = 0.05) 1079> 1080> 1081> ################################################### 1082> ### chunk number 19: monitor1 1083> ################################################### 1084> USIncExp3 <- window(USIncExp2, start = c(1986, 1), end = c(1990,12)) 1085> me.mefp <- monitor(me.mefp) 1086> 1087> 1088> ################################################### 1089> ### chunk number 20: monitor2 1090> ################################################### 1091> USIncExp3 <- window(USIncExp2, start = c(1986, 1)) 1092> me.mefp <- monitor(me.mefp) 1093Break detected at observation # 72 1094> me.mefp 1095Monitoring with ME test (moving estimates test) 1096 1097Initial call: 1098 mefp.formula(formula = ecm.model, type = "ME", data = USIncExp3, alpha = 0.05) 1099 1100Last call: 1101 monitor(obj = me.mefp) 1102 1103Significance level : 0.05 1104Critical value : 3.109524 1105History size : 48 1106Last point evaluated : 182 1107Structural break at : 72 1108 1109Parameter estimate on history : 1110(Intercept) coint.res diff.income 1111 18.9299679 -0.3893141 0.3156597 1112Last parameter estimate : 1113(Intercept) coint.res diff.income 111427.94869106 0.00983451 0.13314662 1115> 1116> 1117> ################################################### 1118> ### chunk number 21: monitor-plot 1119> ################################################### 1120> plot(me.mefp) 1121> 1122> 1123> ################################################### 1124> ### chunk number 22: mefp2 1125> ################################################### 1126> USIncExp3 <- window(USIncExp2, start = c(1986, 1), end = c(1989,12)) 1127> me.efp <- efp(ecm.model, type = "ME", data = USIncExp3, h = 0.5) 1128> me.mefp <- mefp(me.efp, alpha=0.05) 1129> 1130> 1131> ################################################### 1132> ### chunk number 23: monitor3 1133> ################################################### 1134> USIncExp3 <- window(USIncExp2, start = c(1986, 1)) 1135> me.mefp <- monitor(me.mefp) 1136Break detected at observation # 70 1137> 1138> 1139> ################################################### 1140> ### chunk number 24: monitor-plot2 1141> ################################################### 1142> plot(me.mefp) 1143> 1144> 1145> 1146> 1147> cleanEx() 1148> nameEx("boundary.Fstats") 1149> ### * boundary.Fstats 1150> 1151> flush(stderr()); flush(stdout()) 1152> 1153> ### Name: boundary.Fstats 1154> ### Title: Boundary for F Statistics 1155> ### Aliases: boundary.Fstats 1156> ### Keywords: regression 1157> 1158> ### ** Examples 1159> 1160> ## Load dataset "nhtemp" with average yearly temperatures in New Haven 1161> data("nhtemp") 1162> ## plot the data 1163> plot(nhtemp) 1164> 1165> ## test the model null hypothesis that the average temperature remains 1166> ## constant over the years for potential break points between 1941 1167> ## (corresponds to from = 0.5) and 1962 (corresponds to to = 0.85) 1168> ## compute F statistics 1169> fs <- Fstats(nhtemp ~ 1, from = 0.5, to = 0.85) 1170> ## plot the p values without boundary 1171> plot(fs, pval = TRUE, alpha = 0.01) 1172> ## add the boundary in another colour 1173> lines(boundary(fs, pval = TRUE, alpha = 0.01), col = 2) 1174> 1175> 1176> 1177> cleanEx() 1178> nameEx("boundary.efp") 1179> ### * boundary.efp 1180> 1181> flush(stderr()); flush(stdout()) 1182> 1183> ### Name: boundary.efp 1184> ### Title: Boundary for Empirical Fluctuation Processes 1185> ### Aliases: boundary.efp 1186> ### Keywords: regression 1187> 1188> ### ** Examples 1189> 1190> ## Load dataset "nhtemp" with average yearly temperatures in New Haven 1191> data("nhtemp") 1192> ## plot the data 1193> plot(nhtemp) 1194> 1195> ## test the model null hypothesis that the average temperature remains constant 1196> ## over the years 1197> ## compute OLS-CUSUM fluctuation process 1198> temp.cus <- efp(nhtemp ~ 1, type = "OLS-CUSUM") 1199> ## plot the process without boundaries 1200> plot(temp.cus, alpha = 0.01, boundary = FALSE) 1201> ## add the boundaries in another colour 1202> bound <- boundary(temp.cus, alpha = 0.01) 1203> lines(bound, col=4) 1204> lines(-bound, col=4) 1205> 1206> 1207> 1208> cleanEx() 1209> nameEx("boundary.mefp") 1210> ### * boundary.mefp 1211> 1212> flush(stderr()); flush(stdout()) 1213> 1214> ### Name: boundary.mefp 1215> ### Title: Boundary Function for Monitoring of Structural Changes 1216> ### Aliases: boundary.mefp 1217> ### Keywords: regression 1218> 1219> ### ** Examples 1220> 1221> df1 <- data.frame(y=rnorm(300)) 1222> df1[150:300,"y"] <- df1[150:300,"y"]+1 1223> me1 <- mefp(y~1, data=df1[1:50,,drop=FALSE], type="ME", h=1, 1224+ alpha=0.05) 1225> me2 <- monitor(me1, data=df1) 1226Break detected at observation # 183 1227> 1228> plot(me2, boundary=FALSE) 1229> lines(boundary(me2), col="green", lty="44") 1230> 1231> 1232> 1233> cleanEx() 1234> nameEx("breakdates") 1235> ### * breakdates 1236> 1237> flush(stderr()); flush(stdout()) 1238> 1239> ### Name: breakdates 1240> ### Title: Breakdates Corresponding to Breakpoints 1241> ### Aliases: breakdates breakdates.breakpoints 1242> ### breakdates.confint.breakpoints 1243> ### Keywords: regression 1244> 1245> ### ** Examples 1246> 1247> ## Nile data with one breakpoint: the annual flows drop in 1898 1248> ## because the first Ashwan dam was built 1249> data("Nile") 1250> plot(Nile) 1251> 1252> bp.nile <- breakpoints(Nile ~ 1) 1253> summary(bp.nile) 1254 1255 Optimal (m+1)-segment partition: 1256 1257Call: 1258breakpoints.formula(formula = Nile ~ 1) 1259 1260Breakpoints at observation number: 1261 1262m = 1 28 1263m = 2 28 83 1264m = 3 28 68 83 1265m = 4 28 45 68 83 1266m = 5 15 30 45 68 83 1267 1268Corresponding to breakdates: 1269 1270m = 1 1898 1271m = 2 1898 1953 1272m = 3 1898 1938 1953 1273m = 4 1898 1915 1938 1953 1274m = 5 1885 1900 1915 1938 1953 1275 1276Fit: 1277 1278m 0 1 2 3 4 5 1279RSS 2835157 1597457 1552924 1538097 1507888 1659994 1280BIC 1318 1270 1276 1285 1292 1311 1281> plot(bp.nile) 1282> 1283> ## compute breakdates corresponding to the 1284> ## breakpoints of minimum BIC segmentation 1285> breakdates(bp.nile) 1286[1] 1898 1287> 1288> ## confidence intervals 1289> ci.nile <- confint(bp.nile) 1290> breakdates(ci.nile) 1291 2.5 % breakpoints 97.5 % 12921 1895 1898 1902 1293> ci.nile 1294 1295 Confidence intervals for breakpoints 1296 of optimal 2-segment partition: 1297 1298Call: 1299confint.breakpointsfull(object = bp.nile) 1300 1301Breakpoints at observation number: 1302 2.5 % breakpoints 97.5 % 13031 25 28 32 1304 1305Corresponding to breakdates: 1306 2.5 % breakpoints 97.5 % 13071 1895 1898 1902 1308> 1309> plot(Nile) 1310> lines(ci.nile) 1311> 1312> 1313> 1314> cleanEx() 1315> nameEx("breakfactor") 1316> ### * breakfactor 1317> 1318> flush(stderr()); flush(stdout()) 1319> 1320> ### Name: breakfactor 1321> ### Title: Factor Coding of Segmentations 1322> ### Aliases: breakfactor 1323> ### Keywords: regression 1324> 1325> ### ** Examples 1326> 1327> ## Nile data with one breakpoint: the annual flows drop in 1898 1328> ## because the first Ashwan dam was built 1329> data("Nile") 1330> plot(Nile) 1331> 1332> ## compute breakpoints 1333> bp.nile <- breakpoints(Nile ~ 1) 1334> 1335> ## fit and visualize segmented and unsegmented model 1336> fm0 <- lm(Nile ~ 1) 1337> fm1 <- lm(Nile ~ breakfactor(bp.nile, breaks = 1)) 1338> 1339> lines(fitted(fm0), col = 3) 1340> lines(fitted(fm1), col = 4) 1341> lines(bp.nile, breaks = 1) 1342> 1343> 1344> 1345> cleanEx() 1346> nameEx("breakpoints") 1347> ### * breakpoints 1348> 1349> flush(stderr()); flush(stdout()) 1350> 1351> ### Name: breakpoints 1352> ### Title: Dating Breaks 1353> ### Aliases: breakpoints breakpoints.formula breakpoints.breakpointsfull 1354> ### breakpoints.Fstats summary.breakpoints summary.breakpointsfull 1355> ### plot.breakpointsfull plot.summary.breakpointsfull print.breakpoints 1356> ### print.summary.breakpointsfull lines.breakpoints coef.breakpointsfull 1357> ### vcov.breakpointsfull fitted.breakpointsfull residuals.breakpointsfull 1358> ### df.residual.breakpointsfull 1359> ### Keywords: regression 1360> 1361> ### ** Examples 1362> 1363> ## Nile data with one breakpoint: the annual flows drop in 1898 1364> ## because the first Ashwan dam was built 1365> data("Nile") 1366> plot(Nile) 1367> 1368> ## F statistics indicate one breakpoint 1369> fs.nile <- Fstats(Nile ~ 1) 1370> plot(fs.nile) 1371> breakpoints(fs.nile) 1372 1373 Optimal 2-segment partition: 1374 1375Call: 1376breakpoints.Fstats(obj = fs.nile) 1377 1378Breakpoints at observation number: 137928 1380 1381Corresponding to breakdates: 13821898 1383> lines(breakpoints(fs.nile)) 1384> 1385> ## or 1386> bp.nile <- breakpoints(Nile ~ 1) 1387> summary(bp.nile) 1388 1389 Optimal (m+1)-segment partition: 1390 1391Call: 1392breakpoints.formula(formula = Nile ~ 1) 1393 1394Breakpoints at observation number: 1395 1396m = 1 28 1397m = 2 28 83 1398m = 3 28 68 83 1399m = 4 28 45 68 83 1400m = 5 15 30 45 68 83 1401 1402Corresponding to breakdates: 1403 1404m = 1 1898 1405m = 2 1898 1953 1406m = 3 1898 1938 1953 1407m = 4 1898 1915 1938 1953 1408m = 5 1885 1900 1915 1938 1953 1409 1410Fit: 1411 1412m 0 1 2 3 4 5 1413RSS 2835157 1597457 1552924 1538097 1507888 1659994 1414BIC 1318 1270 1276 1285 1292 1311 1415> 1416> ## the BIC also chooses one breakpoint 1417> plot(bp.nile) 1418> breakpoints(bp.nile) 1419 1420 Optimal 2-segment partition: 1421 1422Call: 1423breakpoints.breakpointsfull(obj = bp.nile) 1424 1425Breakpoints at observation number: 142628 1427 1428Corresponding to breakdates: 14291898 1430> 1431> ## fit null hypothesis model and model with 1 breakpoint 1432> fm0 <- lm(Nile ~ 1) 1433> fm1 <- lm(Nile ~ breakfactor(bp.nile, breaks = 1)) 1434> plot(Nile) 1435> lines(ts(fitted(fm0), start = 1871), col = 3) 1436> lines(ts(fitted(fm1), start = 1871), col = 4) 1437> lines(bp.nile) 1438> 1439> ## confidence interval 1440> ci.nile <- confint(bp.nile) 1441> ci.nile 1442 1443 Confidence intervals for breakpoints 1444 of optimal 2-segment partition: 1445 1446Call: 1447confint.breakpointsfull(object = bp.nile) 1448 1449Breakpoints at observation number: 1450 2.5 % breakpoints 97.5 % 14511 25 28 32 1452 1453Corresponding to breakdates: 1454 2.5 % breakpoints 97.5 % 14551 1895 1898 1902 1456> lines(ci.nile) 1457> 1458> 1459> ## UK Seatbelt data: a SARIMA(1,0,0)(1,0,0)_12 model 1460> ## (fitted by OLS) is used and reveals (at least) two 1461> ## breakpoints - one in 1973 associated with the oil crisis and 1462> ## one in 1983 due to the introduction of compulsory 1463> ## wearing of seatbelts in the UK. 1464> data("UKDriverDeaths") 1465> seatbelt <- log10(UKDriverDeaths) 1466> seatbelt <- cbind(seatbelt, lag(seatbelt, k = -1), lag(seatbelt, k = -12)) 1467> colnames(seatbelt) <- c("y", "ylag1", "ylag12") 1468> seatbelt <- window(seatbelt, start = c(1970, 1), end = c(1984,12)) 1469> plot(seatbelt[,"y"], ylab = expression(log[10](casualties))) 1470> 1471> ## testing 1472> re.seat <- efp(y ~ ylag1 + ylag12, data = seatbelt, type = "RE") 1473> plot(re.seat) 1474> 1475> ## dating 1476> bp.seat <- breakpoints(y ~ ylag1 + ylag12, data = seatbelt, h = 0.1) 1477> summary(bp.seat) 1478 1479 Optimal (m+1)-segment partition: 1480 1481Call: 1482breakpoints.formula(formula = y ~ ylag1 + ylag12, h = 0.1, data = seatbelt) 1483 1484Breakpoints at observation number: 1485 1486m = 1 46 1487m = 2 46 157 1488m = 3 46 70 157 1489m = 4 46 70 108 157 1490m = 5 46 70 120 141 160 1491m = 6 46 70 89 108 141 160 1492m = 7 46 70 89 107 125 144 162 1493m = 8 18 46 70 89 107 125 144 162 1494 1495Corresponding to breakdates: 1496 1497m = 1 1973(10) 1498m = 2 1973(10) 1983(1) 1499m = 3 1973(10) 1975(10) 1983(1) 1500m = 4 1973(10) 1975(10) 1978(12) 1983(1) 1501m = 5 1973(10) 1975(10) 1979(12) 1981(9) 1983(4) 1502m = 6 1973(10) 1975(10) 1977(5) 1978(12) 1981(9) 1983(4) 1503m = 7 1973(10) 1975(10) 1977(5) 1978(11) 1980(5) 1981(12) 1983(6) 1504m = 8 1971(6) 1973(10) 1975(10) 1977(5) 1978(11) 1980(5) 1981(12) 1983(6) 1505 1506Fit: 1507 1508m 0 1 2 3 4 5 6 1509RSS 0.3297 0.2967 0.2676 0.2438 0.2395 0.2317 0.2258 1510BIC -602.8611 -601.0539 -598.9042 -594.8774 -577.2905 -562.4880 -546.3632 1511 1512m 7 8 1513RSS 0.2244 0.2231 1514BIC -526.7295 -506.9886 1515> lines(bp.seat, breaks = 2) 1516> 1517> ## minimum BIC partition 1518> plot(bp.seat) 1519> breakpoints(bp.seat) 1520 1521 Optimal 1-segment partition: 1522 1523Call: 1524breakpoints.breakpointsfull(obj = bp.seat) 1525 1526Breakpoints at observation number: 1527NA 1528 1529Corresponding to breakdates: 1530NA 1531> ## the BIC would choose 0 breakpoints although the RE and supF test 1532> ## clearly reject the hypothesis of structural stability. Bai & 1533> ## Perron (2003) report that the BIC has problems in dynamic regressions. 1534> ## due to the shape of the RE process of the F statistics choose two 1535> ## breakpoints and fit corresponding models 1536> bp.seat2 <- breakpoints(bp.seat, breaks = 2) 1537> fm0 <- lm(y ~ ylag1 + ylag12, data = seatbelt) 1538> fm1 <- lm(y ~ breakfactor(bp.seat2)/(ylag1 + ylag12) - 1, data = seatbelt) 1539> 1540> ## plot 1541> plot(seatbelt[,"y"], ylab = expression(log[10](casualties))) 1542> time.seat <- as.vector(time(seatbelt)) 1543> lines(time.seat, fitted(fm0), col = 3) 1544> lines(time.seat, fitted(fm1), col = 4) 1545> lines(bp.seat2) 1546> 1547> ## confidence intervals 1548> ci.seat2 <- confint(bp.seat, breaks = 2) 1549> ci.seat2 1550 1551 Confidence intervals for breakpoints 1552 of optimal 3-segment partition: 1553 1554Call: 1555confint.breakpointsfull(object = bp.seat, breaks = 2) 1556 1557Breakpoints at observation number: 1558 2.5 % breakpoints 97.5 % 15591 33 46 56 15602 144 157 171 1561 1562Corresponding to breakdates: 1563 2.5 % breakpoints 97.5 % 15641 1972(9) 1973(10) 1974(8) 15652 1981(12) 1983(1) 1984(3) 1566> lines(ci.seat2) 1567> 1568> 1569> 1570> cleanEx() 1571> nameEx("catL2BB") 1572> ### * catL2BB 1573> 1574> flush(stderr()); flush(stdout()) 1575> 1576> ### Name: catL2BB 1577> ### Title: Generators for efpFunctionals along Categorical Variables 1578> ### Aliases: catL2BB ordL2BB ordwmax 1579> ### Keywords: regression 1580> 1581> ### ** Examples 1582> 1583> ## artificial data 1584> set.seed(1) 1585> d <- data.frame( 1586+ x = runif(200, -1, 1), 1587+ z = factor(rep(1:4, each = 50)), 1588+ err = rnorm(200) 1589+ ) 1590> d$y <- rep(c(0.5, -0.5), c(150, 50)) * d$x + d$err 1591> 1592> ## empirical fluctuation process 1593> scus <- gefp(y ~ x, data = d, fit = lm, order.by = ~ z) 1594> 1595> ## chi-squared-type test (unordered LM-type test) 1596> LMuo <- catL2BB(scus) 1597> plot(scus, functional = LMuo) 1598> sctest(scus, functional = LMuo) 1599 1600 M-fluctuation test 1601 1602data: scus 1603f(efp) = 12.375, p-value = 0.05411 1604 1605> 1606> ## ordinal maxLM test (with few replications only to save time) 1607> maxLMo <- ordL2BB(scus, nrep = 10000) 1608> plot(scus, functional = maxLMo) 1609> sctest(scus, functional = maxLMo) 1610 1611 M-fluctuation test 1612 1613data: scus 1614f(efp) = 9.0937, p-value = 0.03173 1615 1616> 1617> ## ordinal weighted double maximum test 1618> WDM <- ordwmax(scus) 1619> plot(scus, functional = WDM) 1620> sctest(scus, functional = WDM) 1621 1622 M-fluctuation test 1623 1624data: scus 1625f(efp) = 3.001, p-value = 0.01498 1626 1627> 1628> 1629> 1630> cleanEx() 1631> nameEx("confint.breakpointsfull") 1632> ### * confint.breakpointsfull 1633> 1634> flush(stderr()); flush(stdout()) 1635> 1636> ### Name: confint.breakpointsfull 1637> ### Title: Confidence Intervals for Breakpoints 1638> ### Aliases: confint.breakpointsfull lines.confint.breakpoints 1639> ### print.confint.breakpoints 1640> ### Keywords: regression 1641> 1642> ### ** Examples 1643> 1644> ## Nile data with one breakpoint: the annual flows drop in 1898 1645> ## because the first Ashwan dam was built 1646> data("Nile") 1647> plot(Nile) 1648> 1649> ## dating breaks 1650> bp.nile <- breakpoints(Nile ~ 1) 1651> ci.nile <- confint(bp.nile, breaks = 1) 1652> lines(ci.nile) 1653> 1654> 1655> 1656> cleanEx() 1657> nameEx("durab") 1658> ### * durab 1659> 1660> flush(stderr()); flush(stdout()) 1661> 1662> ### Name: durab 1663> ### Title: US Labor Productivity 1664> ### Aliases: durab 1665> ### Keywords: datasets 1666> 1667> ### ** Examples 1668> 1669> data("durab") 1670> ## use AR(1) model as in Hansen (2001) and Zeileis et al. (2005) 1671> durab.model <- y ~ lag 1672> 1673> ## historical tests 1674> ## OLS-based CUSUM process 1675> ols <- efp(durab.model, data = durab, type = "OLS-CUSUM") 1676> plot(ols) 1677> ## F statistics 1678> fs <- Fstats(durab.model, data = durab, from = 0.1) 1679> plot(fs) 1680> 1681> ## F statistics based on heteroskadisticy-consistent covariance matrix 1682> fsHC <- Fstats(durab.model, data = durab, from = 0.1, 1683+ vcov = function(x, ...) vcovHC(x, type = "HC", ...)) 1684> plot(fsHC) 1685> 1686> ## monitoring 1687> Durab <- window(durab, start=1964, end = c(1979, 12)) 1688> ols.efp <- efp(durab.model, type = "OLS-CUSUM", data = Durab) 1689> newborder <- function(k) 1.723 * k/192 1690> ols.mefp <- mefp(ols.efp, period=2) 1691> ols.mefp2 <- mefp(ols.efp, border=newborder) 1692> Durab <- window(durab, start=1964) 1693> ols.mon <- monitor(ols.mefp) 1694Break detected at observation # 437 1695> ols.mon2 <- monitor(ols.mefp2) 1696Break detected at observation # 416 1697> plot(ols.mon) 1698> lines(boundary(ols.mon2), col = 2) 1699> ## Note: critical value for linear boundary taken from Table III 1700> ## in Zeileis et al. 2005: (1.568 + 1.896)/2 = 1.732 is a linear 1701> ## interpolation between the values for T = 2 and T = 3 at 1702> ## alpha = 0.05. A typo switched 1.732 to 1.723. 1703> 1704> ## dating 1705> bp <- breakpoints(durab.model, data = durab) 1706> summary(bp) 1707 1708 Optimal (m+1)-segment partition: 1709 1710Call: 1711breakpoints.formula(formula = durab.model, data = durab) 1712 1713Breakpoints at observation number: 1714 1715m = 1 418 1716m = 2 221 530 1717m = 3 114 225 530 1718m = 4 114 221 418 531 1719m = 5 114 221 319 418 531 1720 1721Corresponding to breakdates: 1722 1723m = 1 1981(12) 1724m = 2 1965(7) 1991(4) 1725m = 3 1956(8) 1965(11) 1991(4) 1726m = 4 1956(8) 1965(7) 1981(12) 1991(5) 1727m = 5 1956(8) 1965(7) 1973(9) 1981(12) 1991(5) 1728 1729Fit: 1730 1731m 0 1 2 3 4 5 1732RSS 5.586e-02 5.431e-02 5.325e-02 5.220e-02 5.171e-02 5.157e-02 1733BIC -4.221e+03 -4.220e+03 -4.213e+03 -4.207e+03 -4.194e+03 -4.176e+03 1734> plot(summary(bp)) 1735> 1736> plot(ols) 1737> lines(breakpoints(bp, breaks = 1), col = 3) 1738> lines(breakpoints(bp, breaks = 2), col = 4) 1739> plot(fs) 1740> lines(breakpoints(bp, breaks = 1), col = 3) 1741> lines(breakpoints(bp, breaks = 2), col = 4) 1742> 1743> 1744> 1745> cleanEx() 1746> nameEx("efp") 1747> ### * efp 1748> 1749> flush(stderr()); flush(stdout()) 1750> 1751> ### Name: efp 1752> ### Title: Empirical Fluctuation Processes 1753> ### Aliases: efp print.efp 1754> ### Keywords: regression 1755> 1756> ### ** Examples 1757> 1758> ## Nile data with one breakpoint: the annual flows drop in 1898 1759> ## because the first Ashwan dam was built 1760> data("Nile") 1761> plot(Nile) 1762> 1763> ## test the null hypothesis that the annual flow remains constant 1764> ## over the years 1765> ## compute OLS-based CUSUM process and plot 1766> ## with standard and alternative boundaries 1767> ocus.nile <- efp(Nile ~ 1, type = "OLS-CUSUM") 1768> plot(ocus.nile) 1769> plot(ocus.nile, alpha = 0.01, alt.boundary = TRUE) 1770> ## calculate corresponding test statistic 1771> sctest(ocus.nile) 1772 1773 OLS-based CUSUM test 1774 1775data: ocus.nile 1776S0 = 2.9518, p-value = 5.409e-08 1777 1778> 1779> ## UK Seatbelt data: a SARIMA(1,0,0)(1,0,0)_12 model 1780> ## (fitted by OLS) is used and reveals (at least) two 1781> ## breakpoints - one in 1973 associated with the oil crisis and 1782> ## one in 1983 due to the introduction of compulsory 1783> ## wearing of seatbelts in the UK. 1784> data("UKDriverDeaths") 1785> seatbelt <- log10(UKDriverDeaths) 1786> seatbelt <- cbind(seatbelt, lag(seatbelt, k = -1), lag(seatbelt, k = -12)) 1787> colnames(seatbelt) <- c("y", "ylag1", "ylag12") 1788> seatbelt <- window(seatbelt, start = c(1970, 1), end = c(1984,12)) 1789> plot(seatbelt[,"y"], ylab = expression(log[10](casualties))) 1790> 1791> ## use RE process 1792> re.seat <- efp(y ~ ylag1 + ylag12, data = seatbelt, type = "RE") 1793> plot(re.seat) 1794> plot(re.seat, functional = NULL) 1795> sctest(re.seat) 1796 1797 RE test (recursive estimates test) 1798 1799data: re.seat 1800RE = 1.6311, p-value = 0.02904 1801 1802> 1803> 1804> 1805> cleanEx() 1806> nameEx("efpFunctional") 1807> ### * efpFunctional 1808> 1809> flush(stderr()); flush(stdout()) 1810> 1811> ### Name: efpFunctional 1812> ### Title: Functionals for Fluctuation Processes 1813> ### Aliases: efpFunctional simulateBMDist maxBM maxBB maxBMI maxBBI maxL2BB 1814> ### meanL2BB rangeBM rangeBB rangeBMI rangeBBI 1815> ### Keywords: regression 1816> 1817> ### ** Examples 1818> 1819> 1820> data("BostonHomicide") 1821> gcus <- gefp(homicides ~ 1, family = poisson, vcov = kernHAC, 1822+ data = BostonHomicide) 1823> plot(gcus, functional = meanL2BB) 1824> gcus 1825 1826Generalized Empirical M-Fluctuation Process 1827 1828Call: gefp(homicides ~ 1, family = poisson, vcov = kernHAC, data = BostonHomicide) 1829 1830 1831Fitted model: 1832Call: fit(formula = ..1, family = ..2, data = data) 1833 1834Coefficients: 1835(Intercept) 1836 1.017 1837 1838Degrees of Freedom: 76 Total (i.e. Null); 76 Residual 1839Null Deviance: 115.6 1840Residual Deviance: 115.6 AIC: 316.5 1841 1842> sctest(gcus, functional = meanL2BB) 1843 1844 M-fluctuation test 1845 1846data: gcus 1847f(efp) = 0.93375, p-value = 0.005 1848 1849> 1850> y <- rnorm(1000) 1851> x1 <- runif(1000) 1852> x2 <- runif(1000) 1853> 1854> ## supWald statistic computed by Fstats() 1855> fs <- Fstats(y ~ x1 + x2, from = 0.1) 1856> plot(fs) 1857> sctest(fs) 1858 1859 supF test 1860 1861data: fs 1862sup.F = 12.252, p-value = 0.1161 1863 1864> 1865> ## compare with supLM statistic 1866> scus <- gefp(y ~ x1 + x2, fit = lm) 1867> plot(scus, functional = supLM(0.1)) 1868> sctest(scus, functional = supLM(0.1)) 1869 1870 M-fluctuation test 1871 1872data: scus 1873f(efp) = 12.258, p-value = 0.1158 1874 1875> 1876> ## seatbelt data 1877> data("UKDriverDeaths") 1878> seatbelt <- log10(UKDriverDeaths) 1879> seatbelt <- cbind(seatbelt, lag(seatbelt, k = -1), lag(seatbelt, k = -12)) 1880> colnames(seatbelt) <- c("y", "ylag1", "ylag12") 1881> seatbelt <- window(seatbelt, start = c(1970, 1), end = c(1984,12)) 1882> 1883> scus.seat <- gefp(y ~ ylag1 + ylag12, data = seatbelt) 1884> 1885> ## double maximum test 1886> plot(scus.seat) 1887> ## range test 1888> plot(scus.seat, functional = rangeBB) 1889> ## Cramer-von Mises statistic (Nyblom-Hansen test) 1890> plot(scus.seat, functional = meanL2BB) 1891> ## supLM test 1892> plot(scus.seat, functional = supLM(0.1)) 1893> 1894> 1895> 1896> cleanEx() 1897> nameEx("gefp") 1898> ### * gefp 1899> 1900> flush(stderr()); flush(stdout()) 1901> 1902> ### Name: gefp 1903> ### Title: Generalized Empirical M-Fluctuation Processes 1904> ### Aliases: gefp print.gefp sctest.gefp plot.gefp time.gefp print.gefp 1905> ### Keywords: regression 1906> 1907> ### ** Examples 1908> 1909> data("BostonHomicide") 1910> gcus <- gefp(homicides ~ 1, family = poisson, vcov = kernHAC, 1911+ data = BostonHomicide) 1912> plot(gcus, aggregate = FALSE) 1913> gcus 1914 1915Generalized Empirical M-Fluctuation Process 1916 1917Call: gefp(homicides ~ 1, family = poisson, vcov = kernHAC, data = BostonHomicide) 1918 1919 1920Fitted model: 1921Call: fit(formula = ..1, family = ..2, data = data) 1922 1923Coefficients: 1924(Intercept) 1925 1.017 1926 1927Degrees of Freedom: 76 Total (i.e. Null); 76 Residual 1928Null Deviance: 115.6 1929Residual Deviance: 115.6 AIC: 316.5 1930 1931> sctest(gcus) 1932 1933 M-fluctuation test 1934 1935data: gcus 1936f(efp) = 1.669, p-value = 0.007613 1937 1938> 1939> 1940> 1941> cleanEx() 1942> nameEx("logLik.breakpoints") 1943> ### * logLik.breakpoints 1944> 1945> flush(stderr()); flush(stdout()) 1946> 1947> ### Name: logLik.breakpoints 1948> ### Title: Log Likelihood and Information Criteria for Breakpoints 1949> ### Aliases: logLik.breakpoints logLik.breakpointsfull AIC.breakpointsfull 1950> ### Keywords: regression 1951> 1952> ### ** Examples 1953> 1954> ## Nile data with one breakpoint: the annual flows drop in 1898 1955> ## because the first Ashwan dam was built 1956> data("Nile") 1957> plot(Nile) 1958> 1959> bp.nile <- breakpoints(Nile ~ 1) 1960> summary(bp.nile) 1961 1962 Optimal (m+1)-segment partition: 1963 1964Call: 1965breakpoints.formula(formula = Nile ~ 1) 1966 1967Breakpoints at observation number: 1968 1969m = 1 28 1970m = 2 28 83 1971m = 3 28 68 83 1972m = 4 28 45 68 83 1973m = 5 15 30 45 68 83 1974 1975Corresponding to breakdates: 1976 1977m = 1 1898 1978m = 2 1898 1953 1979m = 3 1898 1938 1953 1980m = 4 1898 1915 1938 1953 1981m = 5 1885 1900 1915 1938 1953 1982 1983Fit: 1984 1985m 0 1 2 3 4 5 1986RSS 2835157 1597457 1552924 1538097 1507888 1659994 1987BIC 1318 1270 1276 1285 1292 1311 1988> plot(bp.nile) 1989> 1990> ## BIC of partitions with0 to 5 breakpoints 1991> plot(0:5, AIC(bp.nile, k = log(bp.nile$nobs)), type = "b") 1992> ## AIC 1993> plot(0:5, AIC(bp.nile), type = "b") 1994> 1995> ## BIC, AIC, log likelihood of a single partition 1996> bp.nile1 <- breakpoints(bp.nile, breaks = 1) 1997> AIC(bp.nile1, k = log(bp.nile1$nobs)) 1998[1] 1270.084 1999> AIC(bp.nile1) 2000[1] 1259.663 2001> logLik(bp.nile1) 2002'log Lik.' -625.8315 (df=4) 2003> 2004> 2005> 2006> cleanEx() 2007> nameEx("mefp") 2008> ### * mefp 2009> 2010> flush(stderr()); flush(stdout()) 2011> 2012> ### Name: mefp 2013> ### Title: Monitoring of Empirical Fluctuation Processes 2014> ### Aliases: mefp mefp.formula mefp.efp print.mefp monitor 2015> ### Keywords: regression 2016> 2017> ### ** Examples 2018> 2019> df1 <- data.frame(y=rnorm(300)) 2020> df1[150:300,"y"] <- df1[150:300,"y"]+1 2021> 2022> ## use the first 50 observations as history period 2023> e1 <- efp(y~1, data=df1[1:50,,drop=FALSE], type="ME", h=1) 2024> me1 <- mefp(e1, alpha=0.05) 2025> 2026> ## the same in one function call 2027> me1 <- mefp(y~1, data=df1[1:50,,drop=FALSE], type="ME", h=1, 2028+ alpha=0.05) 2029> 2030> ## monitor the 50 next observations 2031> me2 <- monitor(me1, data=df1[1:100,,drop=FALSE]) 2032> plot(me2) 2033> 2034> # and now monitor on all data 2035> me3 <- monitor(me2, data=df1) 2036Break detected at observation # 183 2037> plot(me3) 2038> 2039> 2040> ## Load dataset "USIncExp" with income and expenditure in the US 2041> ## and choose a suitable subset for the history period 2042> data("USIncExp") 2043> USIncExp3 <- window(USIncExp, start=c(1969,1), end=c(1971,12)) 2044> ## initialize the monitoring with the formula interface 2045> me.mefp <- mefp(expenditure~income, type="ME", rescale=TRUE, 2046+ data=USIncExp3, alpha=0.05) 2047> 2048> ## monitor the new observations for the year 1972 2049> USIncExp3 <- window(USIncExp, start=c(1969,1), end=c(1972,12)) 2050> me.mefp <- monitor(me.mefp) 2051> 2052> ## monitor the new data for the years 1973-1976 2053> USIncExp3 <- window(USIncExp, start=c(1969,1), end=c(1976,12)) 2054> me.mefp <- monitor(me.mefp) 2055Break detected at observation # 58 2056> plot(me.mefp, functional = NULL) 2057> 2058> 2059> 2060> cleanEx() 2061> nameEx("plot.Fstats") 2062> ### * plot.Fstats 2063> 2064> flush(stderr()); flush(stdout()) 2065> 2066> ### Name: plot.Fstats 2067> ### Title: Plot F Statistics 2068> ### Aliases: plot.Fstats lines.Fstats 2069> ### Keywords: hplot 2070> 2071> ### ** Examples 2072> 2073> ## Load dataset "nhtemp" with average yearly temperatures in New Haven 2074> data("nhtemp") 2075> ## plot the data 2076> plot(nhtemp) 2077> 2078> ## test the model null hypothesis that the average temperature remains 2079> ## constant over the years for potential break points between 1941 2080> ## (corresponds to from = 0.5) and 1962 (corresponds to to = 0.85) 2081> ## compute F statistics 2082> fs <- Fstats(nhtemp ~ 1, from = 0.5, to = 0.85) 2083> ## plot the F statistics 2084> plot(fs, alpha = 0.01) 2085> ## and the corresponding p values 2086> plot(fs, pval = TRUE, alpha = 0.01) 2087> ## perform the aveF test 2088> sctest(fs, type = "aveF") 2089 2090 aveF test 2091 2092data: fs 2093ave.F = 10.81, p-value = 2.059e-06 2094 2095> 2096> 2097> 2098> cleanEx() 2099> nameEx("plot.efp") 2100> ### * plot.efp 2101> 2102> flush(stderr()); flush(stdout()) 2103> 2104> ### Name: plot.efp 2105> ### Title: Plot Empirical Fluctuation Process 2106> ### Aliases: plot.efp lines.efp 2107> ### Keywords: hplot 2108> 2109> ### ** Examples 2110> 2111> ## Load dataset "nhtemp" with average yearly temperatures in New Haven 2112> data("nhtemp") 2113> ## plot the data 2114> plot(nhtemp) 2115> 2116> ## test the model null hypothesis that the average temperature remains 2117> ## constant over the years 2118> ## compute Rec-CUSUM fluctuation process 2119> temp.cus <- efp(nhtemp ~ 1) 2120> ## plot the process 2121> plot(temp.cus, alpha = 0.01) 2122> ## and calculate the test statistic 2123> sctest(temp.cus) 2124 2125 Recursive CUSUM test 2126 2127data: temp.cus 2128S = 1.2724, p-value = 0.002902 2129 2130> 2131> ## compute (recursive estimates) fluctuation process 2132> ## with an additional linear trend regressor 2133> lin.trend <- 1:60 2134> temp.me <- efp(nhtemp ~ lin.trend, type = "fluctuation") 2135> ## plot the bivariate process 2136> plot(temp.me, functional = NULL) 2137> ## and perform the corresponding test 2138> sctest(temp.me) 2139 2140 RE test (recursive estimates test) 2141 2142data: temp.me 2143RE = 1.4938, p-value = 0.04558 2144 2145> 2146> 2147> 2148> cleanEx() 2149> nameEx("plot.mefp") 2150> ### * plot.mefp 2151> 2152> flush(stderr()); flush(stdout()) 2153> 2154> ### Name: plot.mefp 2155> ### Title: Plot Methods for mefp Objects 2156> ### Aliases: plot.mefp lines.mefp 2157> ### Keywords: hplot 2158> 2159> ### ** Examples 2160> 2161> df1 <- data.frame(y=rnorm(300)) 2162> df1[150:300,"y"] <- df1[150:300,"y"]+1 2163> me1 <- mefp(y~1, data=df1[1:50,,drop=FALSE], type="ME", h=1, 2164+ alpha=0.05) 2165> me2 <- monitor(me1, data=df1) 2166Break detected at observation # 183 2167> 2168> plot(me2) 2169> 2170> 2171> 2172> cleanEx() 2173> nameEx("recresid") 2174> ### * recresid 2175> 2176> flush(stderr()); flush(stdout()) 2177> 2178> ### Name: recresid 2179> ### Title: Recursive Residuals 2180> ### Aliases: recresid recresid.default recresid.formula recresid.lm 2181> ### Keywords: regression 2182> 2183> ### ** Examples 2184> 2185> x <- rnorm(100) + rep(c(0, 2), each = 50) 2186> rr <- recresid(x ~ 1) 2187> plot(cumsum(rr), type = "l") 2188> 2189> plot(efp(x ~ 1, type = "Rec-CUSUM")) 2190> 2191> 2192> 2193> cleanEx() 2194> nameEx("root.matrix") 2195> ### * root.matrix 2196> 2197> flush(stderr()); flush(stdout()) 2198> 2199> ### Name: root.matrix 2200> ### Title: Root of a Matrix 2201> ### Aliases: root.matrix 2202> ### Keywords: algebra 2203> 2204> ### ** Examples 2205> 2206> X <- matrix(c(1,2,2,8), ncol=2) 2207> test <- root.matrix(X) 2208> ## control results 2209> X 2210 [,1] [,2] 2211[1,] 1 2 2212[2,] 2 8 2213> test %*% test 2214 [,1] [,2] 2215[1,] 1 2 2216[2,] 2 8 2217> 2218> 2219> 2220> cleanEx() 2221> nameEx("scPublications") 2222> ### * scPublications 2223> 2224> flush(stderr()); flush(stdout()) 2225> 2226> ### Name: scPublications 2227> ### Title: Structural Change Publications 2228> ### Aliases: scPublications 2229> ### Keywords: datasets 2230> 2231> ### ** Examples 2232> 2233> ## construct time series: 2234> ## number of sc publications in econometrics/statistics 2235> data("scPublications") 2236> 2237> ## select years from 1987 and 2238> ## `most important' journals 2239> pub <- scPublications 2240> pub <- subset(pub, year > 1986) 2241> tab1 <- table(pub$journal) 2242> nam1 <- names(tab1)[as.vector(tab1) > 9] ## at least 10 papers 2243> tab2 <- sapply(levels(pub$journal), function(x) min(subset(pub, journal == x)$year)) 2244> nam2 <- names(tab2)[as.vector(tab2) < 1991] ## started at least in 1990 2245> nam <- nam1[nam1 %in% nam2] 2246> pub <- subset(pub, as.character(journal) %in% nam) 2247> pub$journal <- factor(pub$journal) 2248> pub_data <- pub 2249> 2250> ## generate time series 2251> pub <- with(pub, tapply(type, year, table)) 2252> pub <- zoo(t(sapply(pub, cbind)), 1987:2006) 2253> colnames(pub) <- levels(pub_data$type) 2254> 2255> ## visualize 2256> plot(pub, ylim = c(0, 35)) 2257> 2258> 2259> 2260> cleanEx() 2261> nameEx("sctest.Fstats") 2262> ### * sctest.Fstats 2263> 2264> flush(stderr()); flush(stdout()) 2265> 2266> ### Name: sctest.Fstats 2267> ### Title: supF-, aveF- and expF-Test 2268> ### Aliases: sctest.Fstats 2269> ### Keywords: htest 2270> 2271> ### ** Examples 2272> 2273> ## Load dataset "nhtemp" with average yearly temperatures in New Haven 2274> data(nhtemp) 2275> ## plot the data 2276> plot(nhtemp) 2277> 2278> ## test the model null hypothesis that the average temperature remains 2279> ## constant over the years for potential break points between 1941 2280> ## (corresponds to from = 0.5) and 1962 (corresponds to to = 0.85) 2281> ## compute F statistics 2282> fs <- Fstats(nhtemp ~ 1, from = 0.5, to = 0.85) 2283> ## plot the F statistics 2284> plot(fs, alpha = 0.01) 2285> ## and the corresponding p values 2286> plot(fs, pval = TRUE, alpha = 0.01) 2287> ## perform the aveF test 2288> sctest(fs, type = "aveF") 2289 2290 aveF test 2291 2292data: fs 2293ave.F = 10.81, p-value = 2.059e-06 2294 2295> 2296> 2297> 2298> cleanEx() 2299> nameEx("sctest.default") 2300> ### * sctest.default 2301> 2302> flush(stderr()); flush(stdout()) 2303> 2304> ### Name: sctest.default 2305> ### Title: Structural Change Tests in Parametric Models 2306> ### Aliases: sctest.default 2307> ### Keywords: htest 2308> 2309> ### ** Examples 2310> 2311> ## Zeileis and Hornik (2007), Section 5.3, Figure 6 2312> data("Grossarl") 2313> m <- glm(cbind(illegitimate, legitimate) ~ 1, family = binomial, data = Grossarl, 2314+ subset = time(fraction) <= 1800) 2315> sctest(m, order.by = 1700:1800, functional = "CvM") 2316 2317 M-fluctuation test 2318 2319data: m 2320f(efp) = 3.5363, p-value = 0.005 2321 2322> 2323> 2324> 2325> cleanEx() 2326> nameEx("sctest.efp") 2327> ### * sctest.efp 2328> 2329> flush(stderr()); flush(stdout()) 2330> 2331> ### Name: sctest.efp 2332> ### Title: Generalized Fluctuation Tests 2333> ### Aliases: sctest.efp 2334> ### Keywords: htest 2335> 2336> ### ** Examples 2337> 2338> ## Load dataset "nhtemp" with average yearly temperatures in New Haven 2339> data("nhtemp") 2340> ## plot the data 2341> plot(nhtemp) 2342> 2343> ## test the model null hypothesis that the average temperature remains 2344> ## constant over the years compute OLS-CUSUM fluctuation process 2345> temp.cus <- efp(nhtemp ~ 1, type = "OLS-CUSUM") 2346> ## plot the process with alternative boundaries 2347> plot(temp.cus, alpha = 0.01, alt.boundary = TRUE) 2348> ## and calculate the test statistic 2349> sctest(temp.cus) 2350 2351 OLS-based CUSUM test 2352 2353data: temp.cus 2354S0 = 2.0728, p-value = 0.0003709 2355 2356> 2357> ## compute moving estimates fluctuation process 2358> temp.me <- efp(nhtemp ~ 1, type = "ME", h = 0.2) 2359> ## plot the process with functional = "max" 2360> plot(temp.me) 2361> ## and perform the corresponding test 2362> sctest(temp.me) 2363 2364 ME test (moving estimates test) 2365 2366data: temp.me 2367ME = 1.5627, p-value = 0.01 2368 2369> 2370> 2371> 2372> cleanEx() 2373> nameEx("sctest.formula") 2374> ### * sctest.formula 2375> 2376> flush(stderr()); flush(stdout()) 2377> 2378> ### Name: sctest.formula 2379> ### Title: Structural Change Tests in Linear Regression Models 2380> ### Aliases: sctest.formula 2381> ### Keywords: htest 2382> 2383> ### ** Examples 2384> 2385> ## Example 7.4 from Greene (1993), "Econometric Analysis" 2386> ## Chow test on Longley data 2387> data("longley") 2388> sctest(Employed ~ Year + GNP.deflator + GNP + Armed.Forces, data = longley, 2389+ type = "Chow", point = 7) 2390 2391 Chow test 2392 2393data: Employed ~ Year + GNP.deflator + GNP + Armed.Forces 2394F = 3.9268, p-value = 0.06307 2395 2396> 2397> ## which is equivalent to segmenting the regression via 2398> fac <- factor(c(rep(1, 7), rep(2, 9))) 2399> fm0 <- lm(Employed ~ Year + GNP.deflator + GNP + Armed.Forces, data = longley) 2400> fm1 <- lm(Employed ~ fac/(Year + GNP.deflator + GNP + Armed.Forces), data = longley) 2401> anova(fm0, fm1) 2402Analysis of Variance Table 2403 2404Model 1: Employed ~ Year + GNP.deflator + GNP + Armed.Forces 2405Model 2: Employed ~ fac/(Year + GNP.deflator + GNP + Armed.Forces) 2406 Res.Df RSS Df Sum of Sq F Pr(>F) 24071 11 4.8987 24082 6 1.1466 5 3.7521 3.9268 0.06307 . 2409--- 2410Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 2411> 2412> ## estimates from Table 7.5 in Greene (1993) 2413> summary(fm0) 2414 2415Call: 2416lm(formula = Employed ~ Year + GNP.deflator + GNP + Armed.Forces, 2417 data = longley) 2418 2419Residuals: 2420 Min 1Q Median 3Q Max 2421-0.9058 -0.3427 -0.1076 0.2168 1.4377 2422 2423Coefficients: 2424 Estimate Std. Error t value Pr(>|t|) 2425(Intercept) 1.169e+03 8.359e+02 1.399 0.18949 2426Year -5.765e-01 4.335e-01 -1.330 0.21049 2427GNP.deflator -1.977e-02 1.389e-01 -0.142 0.88940 2428GNP 6.439e-02 1.995e-02 3.227 0.00805 ** 2429Armed.Forces -1.015e-04 3.086e-03 -0.033 0.97436 2430--- 2431Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 2432 2433Residual standard error: 0.6673 on 11 degrees of freedom 2434Multiple R-squared: 0.9735, Adjusted R-squared: 0.9639 2435F-statistic: 101.1 on 4 and 11 DF, p-value: 1.346e-08 2436 2437> summary(fm1) 2438 2439Call: 2440lm(formula = Employed ~ fac/(Year + GNP.deflator + GNP + Armed.Forces), 2441 data = longley) 2442 2443Residuals: 2444 Min 1Q Median 3Q Max 2445-0.47717 -0.18950 0.02089 0.14836 0.56493 2446 2447Coefficients: 2448 Estimate Std. Error t value Pr(>|t|) 2449(Intercept) 1.678e+03 9.390e+02 1.787 0.12413 2450fac2 2.098e+03 1.786e+03 1.174 0.28473 2451fac1:Year -8.352e-01 4.847e-01 -1.723 0.13563 2452fac2:Year -1.914e+00 7.913e-01 -2.419 0.05194 . 2453fac1:GNP.deflator -1.633e-01 1.762e-01 -0.927 0.38974 2454fac2:GNP.deflator -4.247e-02 2.238e-01 -0.190 0.85576 2455fac1:GNP 9.481e-02 3.815e-02 2.485 0.04747 * 2456fac2:GNP 1.123e-01 2.269e-02 4.951 0.00258 ** 2457fac1:Armed.Forces -2.467e-03 6.965e-03 -0.354 0.73532 2458fac2:Armed.Forces -2.579e-02 1.259e-02 -2.049 0.08635 . 2459--- 2460Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 2461 2462Residual standard error: 0.4372 on 6 degrees of freedom 2463Multiple R-squared: 0.9938, Adjusted R-squared: 0.9845 2464F-statistic: 106.9 on 9 and 6 DF, p-value: 6.28e-06 2465 2466> 2467> 2468> 2469> cleanEx() 2470> nameEx("solveCrossprod") 2471> ### * solveCrossprod 2472> 2473> flush(stderr()); flush(stdout()) 2474> 2475> ### Name: solveCrossprod 2476> ### Title: Inversion of X'X 2477> ### Aliases: solveCrossprod 2478> ### Keywords: algebra 2479> 2480> ### ** Examples 2481> 2482> X <- cbind(1, rnorm(100)) 2483> solveCrossprod(X) 2484 [,1] [,2] 2485[1,] 0.010148448 -0.001363317 2486[2,] -0.001363317 0.012520432 2487> solve(crossprod(X)) 2488 [,1] [,2] 2489[1,] 0.010148448 -0.001363317 2490[2,] -0.001363317 0.012520432 2491> 2492> 2493> 2494> cleanEx() 2495> nameEx("supLM") 2496> ### * supLM 2497> 2498> flush(stderr()); flush(stdout()) 2499> 2500> ### Name: supLM 2501> ### Title: Generators for efpFunctionals along Continuous Variables 2502> ### Aliases: supLM maxMOSUM 2503> ### Keywords: regression 2504> 2505> ### ** Examples 2506> 2507> ## seatbelt data 2508> data("UKDriverDeaths") 2509> seatbelt <- log10(UKDriverDeaths) 2510> seatbelt <- cbind(seatbelt, lag(seatbelt, k = -1), lag(seatbelt, k = -12)) 2511> colnames(seatbelt) <- c("y", "ylag1", "ylag12") 2512> seatbelt <- window(seatbelt, start = c(1970, 1), end = c(1984,12)) 2513> 2514> ## empirical fluctuation process 2515> scus.seat <- gefp(y ~ ylag1 + ylag12, data = seatbelt) 2516> 2517> ## supLM test 2518> plot(scus.seat, functional = supLM(0.1)) 2519> ## MOSUM test 2520> plot(scus.seat, functional = maxMOSUM(0.25)) 2521> ## double maximum test 2522> plot(scus.seat) 2523> ## range test 2524> plot(scus.seat, functional = rangeBB) 2525> ## Cramer-von Mises statistic (Nyblom-Hansen test) 2526> plot(scus.seat, functional = meanL2BB) 2527> 2528> 2529> 2530> ### * <FOOTER> 2531> ### 2532> options(digits = 7L) 2533> base::cat("Time elapsed: ", proc.time() - base::get("ptime", pos = 'CheckExEnv'),"\n") 2534Time elapsed: 43.26 0.124 43.969 0 0 2535> grDevices::dev.off() 2536null device 2537 1 2538> ### 2539> ### Local variables: *** 2540> ### mode: outline-minor *** 2541> ### outline-regexp: "\\(> \\)?### [*]+" *** 2542> ### End: *** 2543> quit('no') 2544