R version 2.9.0 (2009-04-17)
Copyright (C) 2009 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> x <- array(list(337302
+ ,488180
+ ,317756
+ ,318672
+ ,326225
+ ,327532
+ ,338653
+ ,344744
+ ,349420
+ ,520564
+ ,337302
+ ,317756
+ ,318672
+ ,326225
+ ,327532
+ ,338653
+ ,336923
+ ,501492
+ ,349420
+ ,337302
+ ,317756
+ ,318672
+ ,326225
+ ,327532
+ ,330758
+ ,485025
+ ,336923
+ ,349420
+ ,337302
+ ,317756
+ ,318672
+ ,326225
+ ,321002
+ ,464196
+ ,330758
+ ,336923
+ ,349420
+ ,337302
+ ,317756
+ ,318672
+ ,320820
+ ,460170
+ ,321002
+ ,330758
+ ,336923
+ ,349420
+ ,337302
+ ,317756
+ ,327032
+ ,467037
+ ,320820
+ ,321002
+ ,330758
+ ,336923
+ ,349420
+ ,337302
+ ,324047
+ ,460070
+ ,327032
+ ,320820
+ ,321002
+ ,330758
+ ,336923
+ ,349420
+ ,316735
+ ,447988
+ ,324047
+ ,327032
+ ,320820
+ ,321002
+ ,330758
+ ,336923
+ ,315710
+ ,442867
+ ,316735
+ ,324047
+ ,327032
+ ,320820
+ ,321002
+ ,330758
+ ,313427
+ ,436087
+ ,315710
+ ,316735
+ ,324047
+ ,327032
+ ,320820
+ ,321002
+ ,310527
+ ,431328
+ ,313427
+ ,315710
+ ,316735
+ ,324047
+ ,327032
+ ,320820
+ ,330962
+ ,484015
+ ,310527
+ ,313427
+ ,315710
+ ,316735
+ ,324047
+ ,327032
+ ,339015
+ ,509673
+ ,330962
+ ,310527
+ ,313427
+ ,315710
+ ,316735
+ ,324047
+ ,341332
+ ,512927
+ ,339015
+ ,330962
+ ,310527
+ ,313427
+ ,315710
+ ,316735
+ ,339092
+ ,502831
+ ,341332
+ ,339015
+ ,330962
+ ,310527
+ ,313427
+ ,315710
+ ,323308
+ ,470984
+ ,339092
+ ,341332
+ ,339015
+ ,330962
+ ,310527
+ ,313427
+ ,325849
+ ,471067
+ ,323308
+ ,339092
+ ,341332
+ ,339015
+ ,330962
+ ,310527
+ ,330675
+ ,476049
+ ,325849
+ ,323308
+ ,339092
+ ,341332
+ ,339015
+ ,330962
+ ,332225
+ ,474605
+ ,330675
+ ,325849
+ ,323308
+ ,339092
+ ,341332
+ ,339015
+ ,331735
+ ,470439
+ ,332225
+ ,330675
+ ,325849
+ ,323308
+ ,339092
+ ,341332
+ ,328047
+ ,461251
+ ,331735
+ ,332225
+ ,330675
+ ,325849
+ ,323308
+ ,339092
+ ,326165
+ ,454724
+ ,328047
+ ,331735
+ ,332225
+ ,330675
+ ,325849
+ ,323308
+ ,327081
+ ,455626
+ ,326165
+ ,328047
+ ,331735
+ ,332225
+ ,330675
+ ,325849
+ ,346764
+ ,516847
+ ,327081
+ ,326165
+ ,328047
+ ,331735
+ ,332225
+ ,330675
+ ,344190
+ ,525192
+ ,346764
+ ,327081
+ ,326165
+ ,328047
+ ,331735
+ ,332225
+ ,343333
+ ,522975
+ ,344190
+ ,346764
+ ,327081
+ ,326165
+ ,328047
+ ,331735
+ ,345777
+ ,518585
+ ,343333
+ ,344190
+ ,346764
+ ,327081
+ ,326165
+ ,328047
+ ,344094
+ ,509239
+ ,345777
+ ,343333
+ ,344190
+ ,346764
+ ,327081
+ ,326165
+ ,348609
+ ,512238
+ ,344094
+ ,345777
+ ,343333
+ ,344190
+ ,346764
+ ,327081
+ ,354846
+ ,519164
+ ,348609
+ ,344094
+ ,345777
+ ,343333
+ ,344190
+ ,346764
+ ,356427
+ ,517009
+ ,354846
+ ,348609
+ ,344094
+ ,345777
+ ,343333
+ ,344190
+ ,353467
+ ,509933
+ ,356427
+ ,354846
+ ,348609
+ ,344094
+ ,345777
+ ,343333
+ ,355996
+ ,509127
+ ,353467
+ ,356427
+ ,354846
+ ,348609
+ ,344094
+ ,345777
+ ,352487
+ ,500857
+ ,355996
+ ,353467
+ ,356427
+ ,354846
+ ,348609
+ ,344094
+ ,355178
+ ,506971
+ ,352487
+ ,355996
+ ,353467
+ ,356427
+ ,354846
+ ,348609
+ ,374556
+ ,569323
+ ,355178
+ ,352487
+ ,355996
+ ,353467
+ ,356427
+ ,354846
+ ,375021
+ ,579714
+ ,374556
+ ,355178
+ ,352487
+ ,355996
+ ,353467
+ ,356427
+ ,375787
+ ,577992
+ ,375021
+ ,374556
+ ,355178
+ ,352487
+ ,355996
+ ,353467
+ ,372720
+ ,565464
+ ,375787
+ ,375021
+ ,374556
+ ,355178
+ ,352487
+ ,355996
+ ,364431
+ ,547344
+ ,372720
+ ,375787
+ ,375021
+ ,374556
+ ,355178
+ ,352487
+ ,370490
+ ,554788
+ ,364431
+ ,372720
+ ,375787
+ ,375021
+ ,374556
+ ,355178
+ ,376974
+ ,562325
+ ,370490
+ ,364431
+ ,372720
+ ,375787
+ ,375021
+ ,374556
+ ,377632
+ ,560854
+ ,376974
+ ,370490
+ ,364431
+ ,372720
+ ,375787
+ ,375021
+ ,378205
+ ,555332
+ ,377632
+ ,376974
+ ,370490
+ ,364431
+ ,372720
+ ,375787
+ ,370861
+ ,543599
+ ,378205
+ ,377632
+ ,376974
+ ,370490
+ ,364431
+ ,372720
+ ,369167
+ ,536662
+ ,370861
+ ,378205
+ ,377632
+ ,376974
+ ,370490
+ ,364431
+ ,371551
+ ,542722
+ ,369167
+ ,370861
+ ,378205
+ ,377632
+ ,376974
+ ,370490
+ ,382842
+ ,593530
+ ,371551
+ ,369167
+ ,370861
+ ,378205
+ ,377632
+ ,376974
+ ,381903
+ ,610763
+ ,382842
+ ,371551
+ ,369167
+ ,370861
+ ,378205
+ ,377632
+ ,384502
+ ,612613
+ ,381903
+ ,382842
+ ,371551
+ ,369167
+ ,370861
+ ,378205
+ ,392058
+ ,611324
+ ,384502
+ ,381903
+ ,382842
+ ,371551
+ ,369167
+ ,370861
+ ,384359
+ ,594167
+ ,392058
+ ,384502
+ ,381903
+ ,382842
+ ,371551
+ ,369167
+ ,388884
+ ,595454
+ ,384359
+ ,392058
+ ,384502
+ ,381903
+ ,382842
+ ,371551
+ ,386586
+ ,590865
+ ,388884
+ ,384359
+ ,392058
+ ,384502
+ ,381903
+ ,382842
+ ,387495
+ ,589379
+ ,386586
+ ,388884
+ ,384359
+ ,392058
+ ,384502
+ ,381903
+ ,385705
+ ,584428
+ ,387495
+ ,386586
+ ,388884
+ ,384359
+ ,392058
+ ,384502
+ ,378670
+ ,573100
+ ,385705
+ ,387495
+ ,386586
+ ,388884
+ ,384359
+ ,392058
+ ,377367
+ ,567456
+ ,378670
+ ,385705
+ ,387495
+ ,386586
+ ,388884
+ ,384359
+ ,376911
+ ,569028
+ ,377367
+ ,378670
+ ,385705
+ ,387495
+ ,386586
+ ,388884
+ ,389827
+ ,620735
+ ,376911
+ ,377367
+ ,378670
+ ,385705
+ ,387495
+ ,386586
+ ,387820
+ ,628884
+ ,389827
+ ,376911
+ ,377367
+ ,378670
+ ,385705
+ ,387495
+ ,387267
+ ,628232
+ ,387820
+ ,389827
+ ,376911
+ ,377367
+ ,378670
+ ,385705
+ ,380575
+ ,612117
+ ,387267
+ ,387820
+ ,389827
+ ,376911
+ ,377367
+ ,378670
+ ,372402
+ ,595404
+ ,380575
+ ,387267
+ ,387820
+ ,389827
+ ,376911
+ ,377367
+ ,376740
+ ,597141
+ ,372402
+ ,380575
+ ,387267
+ ,387820
+ ,389827
+ ,376911)
+ ,dim=c(8
+ ,66)
+ ,dimnames=list(c('Y'
+ ,'X'
+ ,'yt-1'
+ ,'yt-2'
+ ,'yt-3'
+ ,'yt-4'
+ ,'yt-5'
+ ,'yt-6')
+ ,1:66))
> y <- array(NA,dim=c(8,66),dimnames=list(c('Y','X','yt-1','yt-2','yt-3','yt-4','yt-5','yt-6'),1:66))
> for (i in 1:dim(x)[1])
+ {
+ for (j in 1:dim(x)[2])
+ {
+ y[i,j] <- as.numeric(x[i,j])
+ }
+ }
> par3 = 'Linear Trend'
> par2 = 'Include Monthly Dummies'
> par1 = '1'
> #'GNU S' R Code compiled by R2WASP v. 1.0.44 ()
> #Author: Prof. Dr. P. Wessa
> #To cite this work: AUTHOR(S), (YEAR), YOUR SOFTWARE TITLE (vNUMBER) in Free Statistics Software (v$_version), Office for Research Development and Education, URL http://www.wessa.net/rwasp_YOURPAGE.wasp/
> #Source of accompanying publication: Office for Research, Development, and Education
> #Technical description: Write here your technical program description (don't use hard returns!)
> library(lattice)
> library(lmtest)
Loading required package: zoo
Attaching package: 'zoo'
The following object(s) are masked from package:base :
as.Date.numeric
> n25 <- 25 #minimum number of obs. for Goldfeld-Quandt test
> par1 <- as.numeric(par1)
> x <- t(y)
> k <- length(x[1,])
> n <- length(x[,1])
> x1 <- cbind(x[,par1], x[,1:k!=par1])
> mycolnames <- c(colnames(x)[par1], colnames(x)[1:k!=par1])
> colnames(x1) <- mycolnames #colnames(x)[par1]
> x <- x1
> if (par3 == 'First Differences'){
+ x2 <- array(0, dim=c(n-1,k), dimnames=list(1:(n-1), paste('(1-B)',colnames(x),sep='')))
+ for (i in 1:n-1) {
+ for (j in 1:k) {
+ x2[i,j] <- x[i+1,j] - x[i,j]
+ }
+ }
+ x <- x2
+ }
> if (par2 == 'Include Monthly Dummies'){
+ x2 <- array(0, dim=c(n,11), dimnames=list(1:n, paste('M', seq(1:11), sep ='')))
+ for (i in 1:11){
+ x2[seq(i,n,12),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> if (par2 == 'Include Quarterly Dummies'){
+ x2 <- array(0, dim=c(n,3), dimnames=list(1:n, paste('Q', seq(1:3), sep ='')))
+ for (i in 1:3){
+ x2[seq(i,n,4),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> k <- length(x[1,])
> if (par3 == 'Linear Trend'){
+ x <- cbind(x, c(1:n))
+ colnames(x)[k+1] <- 't'
+ }
> x
Y X yt-1 yt-2 yt-3 yt-4 yt-5 yt-6 M1 M2 M3 M4 M5 M6 M7
1 337302 488180 317756 318672 326225 327532 338653 344744 1 0 0 0 0 0 0
2 349420 520564 337302 317756 318672 326225 327532 338653 0 1 0 0 0 0 0
3 336923 501492 349420 337302 317756 318672 326225 327532 0 0 1 0 0 0 0
4 330758 485025 336923 349420 337302 317756 318672 326225 0 0 0 1 0 0 0
5 321002 464196 330758 336923 349420 337302 317756 318672 0 0 0 0 1 0 0
6 320820 460170 321002 330758 336923 349420 337302 317756 0 0 0 0 0 1 0
7 327032 467037 320820 321002 330758 336923 349420 337302 0 0 0 0 0 0 1
8 324047 460070 327032 320820 321002 330758 336923 349420 0 0 0 0 0 0 0
9 316735 447988 324047 327032 320820 321002 330758 336923 0 0 0 0 0 0 0
10 315710 442867 316735 324047 327032 320820 321002 330758 0 0 0 0 0 0 0
11 313427 436087 315710 316735 324047 327032 320820 321002 0 0 0 0 0 0 0
12 310527 431328 313427 315710 316735 324047 327032 320820 0 0 0 0 0 0 0
13 330962 484015 310527 313427 315710 316735 324047 327032 1 0 0 0 0 0 0
14 339015 509673 330962 310527 313427 315710 316735 324047 0 1 0 0 0 0 0
15 341332 512927 339015 330962 310527 313427 315710 316735 0 0 1 0 0 0 0
16 339092 502831 341332 339015 330962 310527 313427 315710 0 0 0 1 0 0 0
17 323308 470984 339092 341332 339015 330962 310527 313427 0 0 0 0 1 0 0
18 325849 471067 323308 339092 341332 339015 330962 310527 0 0 0 0 0 1 0
19 330675 476049 325849 323308 339092 341332 339015 330962 0 0 0 0 0 0 1
20 332225 474605 330675 325849 323308 339092 341332 339015 0 0 0 0 0 0 0
21 331735 470439 332225 330675 325849 323308 339092 341332 0 0 0 0 0 0 0
22 328047 461251 331735 332225 330675 325849 323308 339092 0 0 0 0 0 0 0
23 326165 454724 328047 331735 332225 330675 325849 323308 0 0 0 0 0 0 0
24 327081 455626 326165 328047 331735 332225 330675 325849 0 0 0 0 0 0 0
25 346764 516847 327081 326165 328047 331735 332225 330675 1 0 0 0 0 0 0
26 344190 525192 346764 327081 326165 328047 331735 332225 0 1 0 0 0 0 0
27 343333 522975 344190 346764 327081 326165 328047 331735 0 0 1 0 0 0 0
28 345777 518585 343333 344190 346764 327081 326165 328047 0 0 0 1 0 0 0
29 344094 509239 345777 343333 344190 346764 327081 326165 0 0 0 0 1 0 0
30 348609 512238 344094 345777 343333 344190 346764 327081 0 0 0 0 0 1 0
31 354846 519164 348609 344094 345777 343333 344190 346764 0 0 0 0 0 0 1
32 356427 517009 354846 348609 344094 345777 343333 344190 0 0 0 0 0 0 0
33 353467 509933 356427 354846 348609 344094 345777 343333 0 0 0 0 0 0 0
34 355996 509127 353467 356427 354846 348609 344094 345777 0 0 0 0 0 0 0
35 352487 500857 355996 353467 356427 354846 348609 344094 0 0 0 0 0 0 0
36 355178 506971 352487 355996 353467 356427 354846 348609 0 0 0 0 0 0 0
37 374556 569323 355178 352487 355996 353467 356427 354846 1 0 0 0 0 0 0
38 375021 579714 374556 355178 352487 355996 353467 356427 0 1 0 0 0 0 0
39 375787 577992 375021 374556 355178 352487 355996 353467 0 0 1 0 0 0 0
40 372720 565464 375787 375021 374556 355178 352487 355996 0 0 0 1 0 0 0
41 364431 547344 372720 375787 375021 374556 355178 352487 0 0 0 0 1 0 0
42 370490 554788 364431 372720 375787 375021 374556 355178 0 0 0 0 0 1 0
43 376974 562325 370490 364431 372720 375787 375021 374556 0 0 0 0 0 0 1
44 377632 560854 376974 370490 364431 372720 375787 375021 0 0 0 0 0 0 0
45 378205 555332 377632 376974 370490 364431 372720 375787 0 0 0 0 0 0 0
46 370861 543599 378205 377632 376974 370490 364431 372720 0 0 0 0 0 0 0
47 369167 536662 370861 378205 377632 376974 370490 364431 0 0 0 0 0 0 0
48 371551 542722 369167 370861 378205 377632 376974 370490 0 0 0 0 0 0 0
49 382842 593530 371551 369167 370861 378205 377632 376974 1 0 0 0 0 0 0
50 381903 610763 382842 371551 369167 370861 378205 377632 0 1 0 0 0 0 0
51 384502 612613 381903 382842 371551 369167 370861 378205 0 0 1 0 0 0 0
52 392058 611324 384502 381903 382842 371551 369167 370861 0 0 0 1 0 0 0
53 384359 594167 392058 384502 381903 382842 371551 369167 0 0 0 0 1 0 0
54 388884 595454 384359 392058 384502 381903 382842 371551 0 0 0 0 0 1 0
55 386586 590865 388884 384359 392058 384502 381903 382842 0 0 0 0 0 0 1
56 387495 589379 386586 388884 384359 392058 384502 381903 0 0 0 0 0 0 0
57 385705 584428 387495 386586 388884 384359 392058 384502 0 0 0 0 0 0 0
58 378670 573100 385705 387495 386586 388884 384359 392058 0 0 0 0 0 0 0
59 377367 567456 378670 385705 387495 386586 388884 384359 0 0 0 0 0 0 0
60 376911 569028 377367 378670 385705 387495 386586 388884 0 0 0 0 0 0 0
61 389827 620735 376911 377367 378670 385705 387495 386586 1 0 0 0 0 0 0
62 387820 628884 389827 376911 377367 378670 385705 387495 0 1 0 0 0 0 0
63 387267 628232 387820 389827 376911 377367 378670 385705 0 0 1 0 0 0 0
64 380575 612117 387267 387820 389827 376911 377367 378670 0 0 0 1 0 0 0
65 372402 595404 380575 387267 387820 389827 376911 377367 0 0 0 0 1 0 0
66 376740 597141 372402 380575 387267 387820 389827 376911 0 0 0 0 0 1 0
M8 M9 M10 M11 t
1 0 0 0 0 1
2 0 0 0 0 2
3 0 0 0 0 3
4 0 0 0 0 4
5 0 0 0 0 5
6 0 0 0 0 6
7 0 0 0 0 7
8 1 0 0 0 8
9 0 1 0 0 9
10 0 0 1 0 10
11 0 0 0 1 11
12 0 0 0 0 12
13 0 0 0 0 13
14 0 0 0 0 14
15 0 0 0 0 15
16 0 0 0 0 16
17 0 0 0 0 17
18 0 0 0 0 18
19 0 0 0 0 19
20 1 0 0 0 20
21 0 1 0 0 21
22 0 0 1 0 22
23 0 0 0 1 23
24 0 0 0 0 24
25 0 0 0 0 25
26 0 0 0 0 26
27 0 0 0 0 27
28 0 0 0 0 28
29 0 0 0 0 29
30 0 0 0 0 30
31 0 0 0 0 31
32 1 0 0 0 32
33 0 1 0 0 33
34 0 0 1 0 34
35 0 0 0 1 35
36 0 0 0 0 36
37 0 0 0 0 37
38 0 0 0 0 38
39 0 0 0 0 39
40 0 0 0 0 40
41 0 0 0 0 41
42 0 0 0 0 42
43 0 0 0 0 43
44 1 0 0 0 44
45 0 1 0 0 45
46 0 0 1 0 46
47 0 0 0 1 47
48 0 0 0 0 48
49 0 0 0 0 49
50 0 0 0 0 50
51 0 0 0 0 51
52 0 0 0 0 52
53 0 0 0 0 53
54 0 0 0 0 54
55 0 0 0 0 55
56 1 0 0 0 56
57 0 1 0 0 57
58 0 0 1 0 58
59 0 0 0 1 59
60 0 0 0 0 60
61 0 0 0 0 61
62 0 0 0 0 62
63 0 0 0 0 63
64 0 0 0 0 64
65 0 0 0 0 65
66 0 0 0 0 66
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X `yt-1` `yt-2` `yt-3` `yt-4`
4.607e+04 6.266e-01 3.366e-01 -5.185e-02 6.985e-02 -6.904e-02
`yt-5` `yt-6` M1 M2 M3 M4
-2.281e-02 -2.434e-01 -1.626e+04 -2.996e+04 -3.010e+04 -2.601e+04
M5 M6 M7 M8 M9 M10
-2.152e+04 -1.486e+04 -9.614e+03 -6.977e+03 -5.770e+03 -3.321e+03
M11 t
-1.455e+03 -4.812e+02
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-4977.5 -1518.9 193.1 1674.9 4247.8
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.607e+04 1.381e+04 3.336 0.001688 **
X 6.266e-01 6.792e-02 9.224 4.98e-12 ***
`yt-1` 3.366e-01 1.138e-01 2.959 0.004864 **
`yt-2` -5.185e-02 1.247e-01 -0.416 0.679455
`yt-3` 6.985e-02 1.247e-01 0.560 0.578199
`yt-4` -6.904e-02 1.233e-01 -0.560 0.578259
`yt-5` -2.281e-02 1.230e-01 -0.186 0.853625
`yt-6` -2.434e-01 8.738e-02 -2.786 0.007733 **
M1 -1.626e+04 3.917e+03 -4.150 0.000142 ***
M2 -2.996e+04 4.203e+03 -7.129 5.84e-09 ***
M3 -3.010e+04 4.208e+03 -7.153 5.37e-09 ***
M4 -2.601e+04 3.692e+03 -7.045 7.80e-09 ***
M5 -2.152e+04 2.664e+03 -8.079 2.26e-10 ***
M6 -1.486e+04 2.857e+03 -5.200 4.48e-06 ***
M7 -9.614e+03 2.293e+03 -4.193 0.000124 ***
M8 -6.977e+03 2.143e+03 -3.256 0.002124 **
M9 -5.770e+03 2.018e+03 -2.859 0.006370 **
M10 -3.321e+03 2.023e+03 -1.642 0.107398
M11 -1.455e+03 1.657e+03 -0.878 0.384687
t -4.812e+02 8.436e+01 -5.704 8.04e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2532 on 46 degrees of freedom
Multiple R-squared: 0.9924, Adjusted R-squared: 0.9893
F-statistic: 316.9 on 19 and 46 DF, p-value: < 2.2e-16
> if (n > n25) {
+ kp3 <- k + 3
+ nmkm3 <- n - k - 3
+ gqarr <- array(NA, dim=c(nmkm3-kp3+1,3))
+ numgqtests <- 0
+ numsignificant1 <- 0
+ numsignificant5 <- 0
+ numsignificant10 <- 0
+ for (mypoint in kp3:nmkm3) {
+ j <- 0
+ numgqtests <- numgqtests + 1
+ for (myalt in c('greater', 'two.sided', 'less')) {
+ j <- j + 1
+ gqarr[mypoint-kp3+1,j] <- gqtest(mylm, point=mypoint, alternative=myalt)$p.value
+ }
+ if (gqarr[mypoint-kp3+1,2] < 0.01) numsignificant1 <- numsignificant1 + 1
+ if (gqarr[mypoint-kp3+1,2] < 0.05) numsignificant5 <- numsignificant5 + 1
+ if (gqarr[mypoint-kp3+1,2] < 0.10) numsignificant10 <- numsignificant10 + 1
+ }
+ gqarr
+ }
[,1] [,2] [,3]
[1,] 0.005865193 0.01173039 0.99413481
[2,] 0.001435930 0.00287186 0.99856407
[3,] 0.807930899 0.38413820 0.19206910
[4,] 0.942260543 0.11547891 0.05773946
[5,] 0.927329161 0.14534168 0.07267084
[6,] 0.879096468 0.24180706 0.12090353
[7,] 0.818528125 0.36294375 0.18147187
[8,] 0.845724894 0.30855021 0.15427511
[9,] 0.817405884 0.36518823 0.18259412
[10,] 0.736671337 0.52665733 0.26332866
[11,] 0.765192860 0.46961428 0.23480714
[12,] 0.705966628 0.58806674 0.29403337
[13,] 0.704689364 0.59062127 0.29531064
[14,] 0.732451760 0.53509648 0.26754824
[15,] 0.856411275 0.28717745 0.14358872
[16,] 0.874780763 0.25043847 0.12521924
[17,] 0.792074806 0.41585039 0.20792519
[18,] 0.691701366 0.61659727 0.30829863
[19,] 0.599700549 0.80059890 0.40029945
[20,] 0.513131438 0.97373712 0.48686856
[21,] 0.363351924 0.72670385 0.63664808
> postscript(file="/var/www/html/rcomp/tmp/10wdn1258562482.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> plot(x[,1], type='l', main='Actuals and Interpolation', ylab='value of Actuals and Interpolation (dots)', xlab='time or index')
> points(x[,1]-mysum$resid)
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/2pzur1258562482.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> plot(mysum$resid, type='b', pch=19, main='Residuals', ylab='value of Residuals', xlab='time or index')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/3lfvm1258562482.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> hist(mysum$resid, main='Residual Histogram', xlab='values of Residuals')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/4lxpp1258562482.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> densityplot(~mysum$resid,col='black',main='Residual Density Plot', xlab='values of Residuals')
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/5mx531258562482.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> qqnorm(mysum$resid, main='Residual Normal Q-Q Plot')
> qqline(mysum$resid)
> grid()
> dev.off()
null device
1
> (myerror <- as.ts(mysum$resid))
Time Series:
Start = 1
End = 66
Frequency = 1
1 2 3 4 5 6
3119.606315 1206.551583 -4977.468062 -1521.700039 -2164.254835 -1109.166039
7 8 9 10 11 12
195.401126 240.250596 -2744.040238 -2391.459010 -3587.503460 -3361.559180
13 14 15 16 17 18
2669.318093 997.676171 -1510.360037 -3329.208319 -2066.687756 -408.511724
19 20 21 22 23 24
334.176667 2101.810760 2470.373194 1749.494870 228.404042 918.112625
25 26 27 28 29 30
9.226773 56.871170 2702.109059 2186.520090 2585.569445 286.469767
31 32 33 34 35 36
316.886676 -1132.466306 -1178.014912 1398.976425 694.810590 1451.296572
37 38 39 40 41 42
-1413.728004 1079.736827 3303.270191 3606.889482 4247.804122 3167.262321
43 44 45 46 47 48
2691.691369 745.420319 3288.131416 198.509642 2489.603516 1921.022195
49 50 51 52 53 54
-627.928718 -2071.434894 581.883261 1958.604316 -920.304396 190.828953
55 56 57 58 59 60
-3538.155839 -1955.015369 -1836.449460 -955.521928 174.685312 -928.872212
61 62 63 64 65 66
-3756.494459 -1269.400856 -99.434412 -2901.105530 -1682.126580 -2126.883278
> postscript(file="/var/www/html/rcomp/tmp/6nm331258562482.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> dum <- cbind(lag(myerror,k=1),myerror)
> dum
Time Series:
Start = 0
End = 66
Frequency = 1
lag(myerror, k = 1) myerror
0 3119.606315 NA
1 1206.551583 3119.606315
2 -4977.468062 1206.551583
3 -1521.700039 -4977.468062
4 -2164.254835 -1521.700039
5 -1109.166039 -2164.254835
6 195.401126 -1109.166039
7 240.250596 195.401126
8 -2744.040238 240.250596
9 -2391.459010 -2744.040238
10 -3587.503460 -2391.459010
11 -3361.559180 -3587.503460
12 2669.318093 -3361.559180
13 997.676171 2669.318093
14 -1510.360037 997.676171
15 -3329.208319 -1510.360037
16 -2066.687756 -3329.208319
17 -408.511724 -2066.687756
18 334.176667 -408.511724
19 2101.810760 334.176667
20 2470.373194 2101.810760
21 1749.494870 2470.373194
22 228.404042 1749.494870
23 918.112625 228.404042
24 9.226773 918.112625
25 56.871170 9.226773
26 2702.109059 56.871170
27 2186.520090 2702.109059
28 2585.569445 2186.520090
29 286.469767 2585.569445
30 316.886676 286.469767
31 -1132.466306 316.886676
32 -1178.014912 -1132.466306
33 1398.976425 -1178.014912
34 694.810590 1398.976425
35 1451.296572 694.810590
36 -1413.728004 1451.296572
37 1079.736827 -1413.728004
38 3303.270191 1079.736827
39 3606.889482 3303.270191
40 4247.804122 3606.889482
41 3167.262321 4247.804122
42 2691.691369 3167.262321
43 745.420319 2691.691369
44 3288.131416 745.420319
45 198.509642 3288.131416
46 2489.603516 198.509642
47 1921.022195 2489.603516
48 -627.928718 1921.022195
49 -2071.434894 -627.928718
50 581.883261 -2071.434894
51 1958.604316 581.883261
52 -920.304396 1958.604316
53 190.828953 -920.304396
54 -3538.155839 190.828953
55 -1955.015369 -3538.155839
56 -1836.449460 -1955.015369
57 -955.521928 -1836.449460
58 174.685312 -955.521928
59 -928.872212 174.685312
60 -3756.494459 -928.872212
61 -1269.400856 -3756.494459
62 -99.434412 -1269.400856
63 -2901.105530 -99.434412
64 -1682.126580 -2901.105530
65 -2126.883278 -1682.126580
66 NA -2126.883278
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 1206.551583 3119.606315
[2,] -4977.468062 1206.551583
[3,] -1521.700039 -4977.468062
[4,] -2164.254835 -1521.700039
[5,] -1109.166039 -2164.254835
[6,] 195.401126 -1109.166039
[7,] 240.250596 195.401126
[8,] -2744.040238 240.250596
[9,] -2391.459010 -2744.040238
[10,] -3587.503460 -2391.459010
[11,] -3361.559180 -3587.503460
[12,] 2669.318093 -3361.559180
[13,] 997.676171 2669.318093
[14,] -1510.360037 997.676171
[15,] -3329.208319 -1510.360037
[16,] -2066.687756 -3329.208319
[17,] -408.511724 -2066.687756
[18,] 334.176667 -408.511724
[19,] 2101.810760 334.176667
[20,] 2470.373194 2101.810760
[21,] 1749.494870 2470.373194
[22,] 228.404042 1749.494870
[23,] 918.112625 228.404042
[24,] 9.226773 918.112625
[25,] 56.871170 9.226773
[26,] 2702.109059 56.871170
[27,] 2186.520090 2702.109059
[28,] 2585.569445 2186.520090
[29,] 286.469767 2585.569445
[30,] 316.886676 286.469767
[31,] -1132.466306 316.886676
[32,] -1178.014912 -1132.466306
[33,] 1398.976425 -1178.014912
[34,] 694.810590 1398.976425
[35,] 1451.296572 694.810590
[36,] -1413.728004 1451.296572
[37,] 1079.736827 -1413.728004
[38,] 3303.270191 1079.736827
[39,] 3606.889482 3303.270191
[40,] 4247.804122 3606.889482
[41,] 3167.262321 4247.804122
[42,] 2691.691369 3167.262321
[43,] 745.420319 2691.691369
[44,] 3288.131416 745.420319
[45,] 198.509642 3288.131416
[46,] 2489.603516 198.509642
[47,] 1921.022195 2489.603516
[48,] -627.928718 1921.022195
[49,] -2071.434894 -627.928718
[50,] 581.883261 -2071.434894
[51,] 1958.604316 581.883261
[52,] -920.304396 1958.604316
[53,] 190.828953 -920.304396
[54,] -3538.155839 190.828953
[55,] -1955.015369 -3538.155839
[56,] -1836.449460 -1955.015369
[57,] -955.521928 -1836.449460
[58,] 174.685312 -955.521928
[59,] -928.872212 174.685312
[60,] -3756.494459 -928.872212
[61,] -1269.400856 -3756.494459
[62,] -99.434412 -1269.400856
[63,] -2901.105530 -99.434412
[64,] -1682.126580 -2901.105530
[65,] -2126.883278 -1682.126580
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 1206.551583 3119.606315
2 -4977.468062 1206.551583
3 -1521.700039 -4977.468062
4 -2164.254835 -1521.700039
5 -1109.166039 -2164.254835
6 195.401126 -1109.166039
7 240.250596 195.401126
8 -2744.040238 240.250596
9 -2391.459010 -2744.040238
10 -3587.503460 -2391.459010
11 -3361.559180 -3587.503460
12 2669.318093 -3361.559180
13 997.676171 2669.318093
14 -1510.360037 997.676171
15 -3329.208319 -1510.360037
16 -2066.687756 -3329.208319
17 -408.511724 -2066.687756
18 334.176667 -408.511724
19 2101.810760 334.176667
20 2470.373194 2101.810760
21 1749.494870 2470.373194
22 228.404042 1749.494870
23 918.112625 228.404042
24 9.226773 918.112625
25 56.871170 9.226773
26 2702.109059 56.871170
27 2186.520090 2702.109059
28 2585.569445 2186.520090
29 286.469767 2585.569445
30 316.886676 286.469767
31 -1132.466306 316.886676
32 -1178.014912 -1132.466306
33 1398.976425 -1178.014912
34 694.810590 1398.976425
35 1451.296572 694.810590
36 -1413.728004 1451.296572
37 1079.736827 -1413.728004
38 3303.270191 1079.736827
39 3606.889482 3303.270191
40 4247.804122 3606.889482
41 3167.262321 4247.804122
42 2691.691369 3167.262321
43 745.420319 2691.691369
44 3288.131416 745.420319
45 198.509642 3288.131416
46 2489.603516 198.509642
47 1921.022195 2489.603516
48 -627.928718 1921.022195
49 -2071.434894 -627.928718
50 581.883261 -2071.434894
51 1958.604316 581.883261
52 -920.304396 1958.604316
53 190.828953 -920.304396
54 -3538.155839 190.828953
55 -1955.015369 -3538.155839
56 -1836.449460 -1955.015369
57 -955.521928 -1836.449460
58 174.685312 -955.521928
59 -928.872212 174.685312
60 -3756.494459 -928.872212
61 -1269.400856 -3756.494459
62 -99.434412 -1269.400856
63 -2901.105530 -99.434412
64 -1682.126580 -2901.105530
65 -2126.883278 -1682.126580
> plot(z,main=paste('Residual Lag plot, lowess, and regression line'), ylab='values of Residuals', xlab='lagged values of Residuals')
> lines(lowess(z))
> abline(lm(z))
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/71mgr1258562482.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> acf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Autocorrelation Function')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/86l3l1258562482.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> pacf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Partial Autocorrelation Function')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/9jlrt1258562482.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> opar <- par(mfrow = c(2,2), oma = c(0, 0, 1.1, 0))
> plot(mylm, las = 1, sub='Residual Diagnostics')
> par(opar)
> dev.off()
null device
1
> if (n > n25) {
+ postscript(file="/var/www/html/rcomp/tmp/104r0l1258562482.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
+ plot(kp3:nmkm3,gqarr[,2], main='Goldfeld-Quandt test',ylab='2-sided p-value',xlab='breakpoint')
+ grid()
+ dev.off()
+ }
null device
1
>
> #Note: the /var/www/html/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/html/rcomp/createtable")
>
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Estimated Regression Equation', 1, TRUE)
> a<-table.row.end(a)
> myeq <- colnames(x)[1]
> myeq <- paste(myeq, '[t] = ', sep='')
> for (i in 1:k){
+ if (mysum$coefficients[i,1] > 0) myeq <- paste(myeq, '+', '')
+ myeq <- paste(myeq, mysum$coefficients[i,1], sep=' ')
+ if (rownames(mysum$coefficients)[i] != '(Intercept)') {
+ myeq <- paste(myeq, rownames(mysum$coefficients)[i], sep='')
+ if (rownames(mysum$coefficients)[i] != 't') myeq <- paste(myeq, '[t]', sep='')
+ }
+ }
> myeq <- paste(myeq, ' + e[t]')
> a<-table.row.start(a)
> a<-table.element(a, myeq)
> a<-table.row.end(a)
> a<-table.end(a)
> table.save(a,file="/var/www/html/rcomp/tmp/11zo341258562482.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a,hyperlink('http://www.xycoon.com/ols1.htm','Multiple Linear Regression - Ordinary Least Squares',''), 6, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a,'Variable',header=TRUE)
> a<-table.element(a,'Parameter',header=TRUE)
> a<-table.element(a,'S.D.',header=TRUE)
> a<-table.element(a,'T-STAT
H0: parameter = 0',header=TRUE)
> a<-table.element(a,'2-tail p-value',header=TRUE)
> a<-table.element(a,'1-tail p-value',header=TRUE)
> a<-table.row.end(a)
> for (i in 1:k){
+ a<-table.row.start(a)
+ a<-table.element(a,rownames(mysum$coefficients)[i],header=TRUE)
+ a<-table.element(a,mysum$coefficients[i,1])
+ a<-table.element(a, round(mysum$coefficients[i,2],6))
+ a<-table.element(a, round(mysum$coefficients[i,3],4))
+ a<-table.element(a, round(mysum$coefficients[i,4],6))
+ a<-table.element(a, round(mysum$coefficients[i,4]/2,6))
+ a<-table.row.end(a)
+ }
> a<-table.end(a)
> table.save(a,file="/var/www/html/rcomp/tmp/12j0gj1258562482.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Regression Statistics', 2, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple R',1,TRUE)
> a<-table.element(a, sqrt(mysum$r.squared))
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'R-squared',1,TRUE)
> a<-table.element(a, mysum$r.squared)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Adjusted R-squared',1,TRUE)
> a<-table.element(a, mysum$adj.r.squared)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (value)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[1])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (DF numerator)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[2])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (DF denominator)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[3])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'p-value',1,TRUE)
> a<-table.element(a, 1-pf(mysum$fstatistic[1],mysum$fstatistic[2],mysum$fstatistic[3]))
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Residual Statistics', 2, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Residual Standard Deviation',1,TRUE)
> a<-table.element(a, mysum$sigma)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Sum Squared Residuals',1,TRUE)
> a<-table.element(a, sum(myerror*myerror))
> a<-table.row.end(a)
> a<-table.end(a)
> table.save(a,file="/var/www/html/rcomp/tmp/139onv1258562482.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Actuals, Interpolation, and Residuals', 4, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Time or Index', 1, TRUE)
> a<-table.element(a, 'Actuals', 1, TRUE)
> a<-table.element(a, 'Interpolation
Forecast', 1, TRUE)
> a<-table.element(a, 'Residuals
Prediction Error', 1, TRUE)
> a<-table.row.end(a)
> for (i in 1:n) {
+ a<-table.row.start(a)
+ a<-table.element(a,i, 1, TRUE)
+ a<-table.element(a,x[i])
+ a<-table.element(a,x[i]-mysum$resid[i])
+ a<-table.element(a,mysum$resid[i])
+ a<-table.row.end(a)
+ }
> a<-table.end(a)
> table.save(a,file="/var/www/html/rcomp/tmp/14m8yg1258562482.tab")
> if (n > n25) {
+ a<-table.start()
+ a<-table.row.start(a)
+ a<-table.element(a,'Goldfeld-Quandt test for Heteroskedasticity',4,TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'p-values',header=TRUE)
+ a<-table.element(a,'Alternative Hypothesis',3,header=TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'breakpoint index',header=TRUE)
+ a<-table.element(a,'greater',header=TRUE)
+ a<-table.element(a,'2-sided',header=TRUE)
+ a<-table.element(a,'less',header=TRUE)
+ a<-table.row.end(a)
+ for (mypoint in kp3:nmkm3) {
+ a<-table.row.start(a)
+ a<-table.element(a,mypoint,header=TRUE)
+ a<-table.element(a,gqarr[mypoint-kp3+1,1])
+ a<-table.element(a,gqarr[mypoint-kp3+1,2])
+ a<-table.element(a,gqarr[mypoint-kp3+1,3])
+ a<-table.row.end(a)
+ }
+ a<-table.end(a)
+ table.save(a,file="/var/www/html/rcomp/tmp/15uj851258562482.tab")
+ a<-table.start()
+ a<-table.row.start(a)
+ a<-table.element(a,'Meta Analysis of Goldfeld-Quandt test for Heteroskedasticity',4,TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'Description',header=TRUE)
+ a<-table.element(a,'# significant tests',header=TRUE)
+ a<-table.element(a,'% significant tests',header=TRUE)
+ a<-table.element(a,'OK/NOK',header=TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'1% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant1)
+ a<-table.element(a,numsignificant1/numgqtests)
+ if (numsignificant1/numgqtests < 0.01) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'5% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant5)
+ a<-table.element(a,numsignificant5/numgqtests)
+ if (numsignificant5/numgqtests < 0.05) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'10% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant10)
+ a<-table.element(a,numsignificant10/numgqtests)
+ if (numsignificant10/numgqtests < 0.1) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.end(a)
+ table.save(a,file="/var/www/html/rcomp/tmp/1621a71258562482.tab")
+ }
>
> system("convert tmp/10wdn1258562482.ps tmp/10wdn1258562482.png")
> system("convert tmp/2pzur1258562482.ps tmp/2pzur1258562482.png")
> system("convert tmp/3lfvm1258562482.ps tmp/3lfvm1258562482.png")
> system("convert tmp/4lxpp1258562482.ps tmp/4lxpp1258562482.png")
> system("convert tmp/5mx531258562482.ps tmp/5mx531258562482.png")
> system("convert tmp/6nm331258562482.ps tmp/6nm331258562482.png")
> system("convert tmp/71mgr1258562482.ps tmp/71mgr1258562482.png")
> system("convert tmp/86l3l1258562482.ps tmp/86l3l1258562482.png")
> system("convert tmp/9jlrt1258562482.ps tmp/9jlrt1258562482.png")
> system("convert tmp/104r0l1258562482.ps tmp/104r0l1258562482.png")
>
>
> proc.time()
user system elapsed
2.466 1.546 2.840