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(600969
+ ,586840
+ ,671833
+ ,654294
+ ,631923
+ ,625568
+ ,600969
+ ,586840
+ ,671833
+ ,654294
+ ,558110
+ ,625568
+ ,600969
+ ,586840
+ ,671833
+ ,630577
+ ,558110
+ ,625568
+ ,600969
+ ,586840
+ ,628654
+ ,630577
+ ,558110
+ ,625568
+ ,600969
+ ,603184
+ ,628654
+ ,630577
+ ,558110
+ ,625568
+ ,656255
+ ,603184
+ ,628654
+ ,630577
+ ,558110
+ ,600730
+ ,656255
+ ,603184
+ ,628654
+ ,630577
+ ,670326
+ ,600730
+ ,656255
+ ,603184
+ ,628654
+ ,678423
+ ,670326
+ ,600730
+ ,656255
+ ,603184
+ ,641502
+ ,678423
+ ,670326
+ ,600730
+ ,656255
+ ,625311
+ ,641502
+ ,678423
+ ,670326
+ ,600730
+ ,628177
+ ,625311
+ ,641502
+ ,678423
+ ,670326
+ ,589767
+ ,628177
+ ,625311
+ ,641502
+ ,678423
+ ,582471
+ ,589767
+ ,628177
+ ,625311
+ ,641502
+ ,636248
+ ,582471
+ ,589767
+ ,628177
+ ,625311
+ ,599885
+ ,636248
+ ,582471
+ ,589767
+ ,628177
+ ,621694
+ ,599885
+ ,636248
+ ,582471
+ ,589767
+ ,637406
+ ,621694
+ ,599885
+ ,636248
+ ,582471
+ ,595994
+ ,637406
+ ,621694
+ ,599885
+ ,636248
+ ,696308
+ ,595994
+ ,637406
+ ,621694
+ ,599885
+ ,674201
+ ,696308
+ ,595994
+ ,637406
+ ,621694
+ ,648861
+ ,674201
+ ,696308
+ ,595994
+ ,637406
+ ,649605
+ ,648861
+ ,674201
+ ,696308
+ ,595994
+ ,672392
+ ,649605
+ ,648861
+ ,674201
+ ,696308
+ ,598396
+ ,672392
+ ,649605
+ ,648861
+ ,674201
+ ,613177
+ ,598396
+ ,672392
+ ,649605
+ ,648861
+ ,638104
+ ,613177
+ ,598396
+ ,672392
+ ,649605
+ ,615632
+ ,638104
+ ,613177
+ ,598396
+ ,672392
+ ,634465
+ ,615632
+ ,638104
+ ,613177
+ ,598396
+ ,638686
+ ,634465
+ ,615632
+ ,638104
+ ,613177
+ ,604243
+ ,638686
+ ,634465
+ ,615632
+ ,638104
+ ,706669
+ ,604243
+ ,638686
+ ,634465
+ ,615632
+ ,677185
+ ,706669
+ ,604243
+ ,638686
+ ,634465
+ ,644328
+ ,677185
+ ,706669
+ ,604243
+ ,638686
+ ,644825
+ ,644328
+ ,677185
+ ,706669
+ ,604243
+ ,605707
+ ,644825
+ ,644328
+ ,677185
+ ,706669
+ ,600136
+ ,605707
+ ,644825
+ ,644328
+ ,677185
+ ,612166
+ ,600136
+ ,605707
+ ,644825
+ ,644328
+ ,599659
+ ,612166
+ ,600136
+ ,605707
+ ,644825
+ ,634210
+ ,599659
+ ,612166
+ ,600136
+ ,605707
+ ,618234
+ ,634210
+ ,599659
+ ,612166
+ ,600136
+ ,613576
+ ,618234
+ ,634210
+ ,599659
+ ,612166
+ ,627200
+ ,613576
+ ,618234
+ ,634210
+ ,599659
+ ,668973
+ ,627200
+ ,613576
+ ,618234
+ ,634210
+ ,651479
+ ,668973
+ ,627200
+ ,613576
+ ,618234
+ ,619661
+ ,651479
+ ,668973
+ ,627200
+ ,613576
+ ,644260
+ ,619661
+ ,651479
+ ,668973
+ ,627200
+ ,579936
+ ,644260
+ ,619661
+ ,651479
+ ,668973
+ ,601752
+ ,579936
+ ,644260
+ ,619661
+ ,651479
+ ,595376
+ ,601752
+ ,579936
+ ,644260
+ ,619661
+ ,588902
+ ,595376
+ ,601752
+ ,579936
+ ,644260
+ ,634341
+ ,588902
+ ,595376
+ ,601752
+ ,579936
+ ,594305
+ ,634341
+ ,588902
+ ,595376
+ ,601752
+ ,606200
+ ,594305
+ ,634341
+ ,588902
+ ,595376
+ ,610926
+ ,606200
+ ,594305
+ ,634341
+ ,588902
+ ,633685
+ ,610926
+ ,606200
+ ,594305
+ ,634341
+ ,639696
+ ,633685
+ ,610926
+ ,606200
+ ,594305
+ ,659451
+ ,639696
+ ,633685
+ ,610926
+ ,606200
+ ,593248
+ ,659451
+ ,639696
+ ,633685
+ ,610926
+ ,606677
+ ,593248
+ ,659451
+ ,639696
+ ,633685
+ ,599434
+ ,606677
+ ,593248
+ ,659451
+ ,639696
+ ,569578
+ ,599434
+ ,606677
+ ,593248
+ ,659451
+ ,629873
+ ,569578
+ ,599434
+ ,606677
+ ,593248
+ ,613438
+ ,629873
+ ,569578
+ ,599434
+ ,606677
+ ,604172
+ ,613438
+ ,629873
+ ,569578
+ ,599434
+ ,658328
+ ,604172
+ ,613438
+ ,629873
+ ,569578
+ ,612633
+ ,658328
+ ,604172
+ ,613438
+ ,629873
+ ,707372
+ ,612633
+ ,658328
+ ,604172
+ ,613438
+ ,739770
+ ,707372
+ ,612633
+ ,658328
+ ,604172
+ ,777535
+ ,739770
+ ,707372
+ ,612633
+ ,658328
+ ,685030
+ ,777535
+ ,739770
+ ,707372
+ ,612633
+ ,730234
+ ,685030
+ ,777535
+ ,739770
+ ,707372
+ ,714154
+ ,730234
+ ,685030
+ ,777535
+ ,739770
+ ,630872
+ ,714154
+ ,730234
+ ,685030
+ ,777535
+ ,719492
+ ,630872
+ ,714154
+ ,730234
+ ,685030
+ ,677023
+ ,719492
+ ,630872
+ ,714154
+ ,730234
+ ,679272
+ ,677023
+ ,719492
+ ,630872
+ ,714154
+ ,718317
+ ,679272
+ ,677023
+ ,719492
+ ,630872
+ ,645672
+ ,718317
+ ,679272
+ ,677023
+ ,719492)
+ ,dim=c(5
+ ,80)
+ ,dimnames=list(c('yt'
+ ,'yt-1'
+ ,'yt-2'
+ ,'yt-3'
+ ,'yt-4')
+ ,1:80))
> y <- array(NA,dim=c(5,80),dimnames=list(c('yt','yt-1','yt-2','yt-3','yt-4'),1:80))
> 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
yt yt-1 yt-2 yt-3 yt-4 M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 600969 586840 671833 654294 631923 1 0 0 0 0 0 0 0 0 0 0 1
2 625568 600969 586840 671833 654294 0 1 0 0 0 0 0 0 0 0 0 2
3 558110 625568 600969 586840 671833 0 0 1 0 0 0 0 0 0 0 0 3
4 630577 558110 625568 600969 586840 0 0 0 1 0 0 0 0 0 0 0 4
5 628654 630577 558110 625568 600969 0 0 0 0 1 0 0 0 0 0 0 5
6 603184 628654 630577 558110 625568 0 0 0 0 0 1 0 0 0 0 0 6
7 656255 603184 628654 630577 558110 0 0 0 0 0 0 1 0 0 0 0 7
8 600730 656255 603184 628654 630577 0 0 0 0 0 0 0 1 0 0 0 8
9 670326 600730 656255 603184 628654 0 0 0 0 0 0 0 0 1 0 0 9
10 678423 670326 600730 656255 603184 0 0 0 0 0 0 0 0 0 1 0 10
11 641502 678423 670326 600730 656255 0 0 0 0 0 0 0 0 0 0 1 11
12 625311 641502 678423 670326 600730 0 0 0 0 0 0 0 0 0 0 0 12
13 628177 625311 641502 678423 670326 1 0 0 0 0 0 0 0 0 0 0 13
14 589767 628177 625311 641502 678423 0 1 0 0 0 0 0 0 0 0 0 14
15 582471 589767 628177 625311 641502 0 0 1 0 0 0 0 0 0 0 0 15
16 636248 582471 589767 628177 625311 0 0 0 1 0 0 0 0 0 0 0 16
17 599885 636248 582471 589767 628177 0 0 0 0 1 0 0 0 0 0 0 17
18 621694 599885 636248 582471 589767 0 0 0 0 0 1 0 0 0 0 0 18
19 637406 621694 599885 636248 582471 0 0 0 0 0 0 1 0 0 0 0 19
20 595994 637406 621694 599885 636248 0 0 0 0 0 0 0 1 0 0 0 20
21 696308 595994 637406 621694 599885 0 0 0 0 0 0 0 0 1 0 0 21
22 674201 696308 595994 637406 621694 0 0 0 0 0 0 0 0 0 1 0 22
23 648861 674201 696308 595994 637406 0 0 0 0 0 0 0 0 0 0 1 23
24 649605 648861 674201 696308 595994 0 0 0 0 0 0 0 0 0 0 0 24
25 672392 649605 648861 674201 696308 1 0 0 0 0 0 0 0 0 0 0 25
26 598396 672392 649605 648861 674201 0 1 0 0 0 0 0 0 0 0 0 26
27 613177 598396 672392 649605 648861 0 0 1 0 0 0 0 0 0 0 0 27
28 638104 613177 598396 672392 649605 0 0 0 1 0 0 0 0 0 0 0 28
29 615632 638104 613177 598396 672392 0 0 0 0 1 0 0 0 0 0 0 29
30 634465 615632 638104 613177 598396 0 0 0 0 0 1 0 0 0 0 0 30
31 638686 634465 615632 638104 613177 0 0 0 0 0 0 1 0 0 0 0 31
32 604243 638686 634465 615632 638104 0 0 0 0 0 0 0 1 0 0 0 32
33 706669 604243 638686 634465 615632 0 0 0 0 0 0 0 0 1 0 0 33
34 677185 706669 604243 638686 634465 0 0 0 0 0 0 0 0 0 1 0 34
35 644328 677185 706669 604243 638686 0 0 0 0 0 0 0 0 0 0 1 35
36 644825 644328 677185 706669 604243 0 0 0 0 0 0 0 0 0 0 0 36
37 605707 644825 644328 677185 706669 1 0 0 0 0 0 0 0 0 0 0 37
38 600136 605707 644825 644328 677185 0 1 0 0 0 0 0 0 0 0 0 38
39 612166 600136 605707 644825 644328 0 0 1 0 0 0 0 0 0 0 0 39
40 599659 612166 600136 605707 644825 0 0 0 1 0 0 0 0 0 0 0 40
41 634210 599659 612166 600136 605707 0 0 0 0 1 0 0 0 0 0 0 41
42 618234 634210 599659 612166 600136 0 0 0 0 0 1 0 0 0 0 0 42
43 613576 618234 634210 599659 612166 0 0 0 0 0 0 1 0 0 0 0 43
44 627200 613576 618234 634210 599659 0 0 0 0 0 0 0 1 0 0 0 44
45 668973 627200 613576 618234 634210 0 0 0 0 0 0 0 0 1 0 0 45
46 651479 668973 627200 613576 618234 0 0 0 0 0 0 0 0 0 1 0 46
47 619661 651479 668973 627200 613576 0 0 0 0 0 0 0 0 0 0 1 47
48 644260 619661 651479 668973 627200 0 0 0 0 0 0 0 0 0 0 0 48
49 579936 644260 619661 651479 668973 1 0 0 0 0 0 0 0 0 0 0 49
50 601752 579936 644260 619661 651479 0 1 0 0 0 0 0 0 0 0 0 50
51 595376 601752 579936 644260 619661 0 0 1 0 0 0 0 0 0 0 0 51
52 588902 595376 601752 579936 644260 0 0 0 1 0 0 0 0 0 0 0 52
53 634341 588902 595376 601752 579936 0 0 0 0 1 0 0 0 0 0 0 53
54 594305 634341 588902 595376 601752 0 0 0 0 0 1 0 0 0 0 0 54
55 606200 594305 634341 588902 595376 0 0 0 0 0 0 1 0 0 0 0 55
56 610926 606200 594305 634341 588902 0 0 0 0 0 0 0 1 0 0 0 56
57 633685 610926 606200 594305 634341 0 0 0 0 0 0 0 0 1 0 0 57
58 639696 633685 610926 606200 594305 0 0 0 0 0 0 0 0 0 1 0 58
59 659451 639696 633685 610926 606200 0 0 0 0 0 0 0 0 0 0 1 59
60 593248 659451 639696 633685 610926 0 0 0 0 0 0 0 0 0 0 0 60
61 606677 593248 659451 639696 633685 1 0 0 0 0 0 0 0 0 0 0 61
62 599434 606677 593248 659451 639696 0 1 0 0 0 0 0 0 0 0 0 62
63 569578 599434 606677 593248 659451 0 0 1 0 0 0 0 0 0 0 0 63
64 629873 569578 599434 606677 593248 0 0 0 1 0 0 0 0 0 0 0 64
65 613438 629873 569578 599434 606677 0 0 0 0 1 0 0 0 0 0 0 65
66 604172 613438 629873 569578 599434 0 0 0 0 0 1 0 0 0 0 0 66
67 658328 604172 613438 629873 569578 0 0 0 0 0 0 1 0 0 0 0 67
68 612633 658328 604172 613438 629873 0 0 0 0 0 0 0 1 0 0 0 68
69 707372 612633 658328 604172 613438 0 0 0 0 0 0 0 0 1 0 0 69
70 739770 707372 612633 658328 604172 0 0 0 0 0 0 0 0 0 1 0 70
71 777535 739770 707372 612633 658328 0 0 0 0 0 0 0 0 0 0 1 71
72 685030 777535 739770 707372 612633 0 0 0 0 0 0 0 0 0 0 0 72
73 730234 685030 777535 739770 707372 1 0 0 0 0 0 0 0 0 0 0 73
74 714154 730234 685030 777535 739770 0 1 0 0 0 0 0 0 0 0 0 74
75 630872 714154 730234 685030 777535 0 0 1 0 0 0 0 0 0 0 0 75
76 719492 630872 714154 730234 685030 0 0 0 1 0 0 0 0 0 0 0 76
77 677023 719492 630872 714154 730234 0 0 0 0 1 0 0 0 0 0 0 77
78 679272 677023 719492 630872 714154 0 0 0 0 0 1 0 0 0 0 0 78
79 718317 679272 677023 719492 630872 0 0 0 0 0 0 1 0 0 0 0 79
80 645672 718317 679272 677023 719492 0 0 0 0 0 0 0 1 0 0 0 80
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) `yt-1` `yt-2` `yt-3` `yt-4` M1
2.812e+04 2.242e-01 3.652e-01 6.027e-01 -3.279e-01 2.894e+04
M2 M3 M4 M5 M6 M7
3.209e+04 2.907e+04 6.888e+04 7.075e+04 6.037e+04 5.658e+04
M8 M9 M10 M11 t
3.692e+04 1.135e+05 8.686e+04 7.303e+04 1.263e+02
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-54659 -8971 -2606 7418 89814
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.812e+04 6.043e+04 0.465 0.643318
`yt-1` 2.243e-01 1.199e-01 1.870 0.066126 .
`yt-2` 3.652e-01 9.354e-02 3.904 0.000233 ***
`yt-3` 6.027e-01 9.570e-02 6.298 3.29e-08 ***
`yt-4` -3.279e-01 1.198e-01 -2.736 0.008076 **
M1 2.894e+04 1.567e+04 1.848 0.069366 .
M2 3.209e+04 1.653e+04 1.941 0.056687 .
M3 2.907e+04 1.671e+04 1.739 0.086896 .
M4 6.888e+04 1.593e+04 4.323 5.59e-05 ***
M5 7.075e+04 1.510e+04 4.685 1.54e-05 ***
M6 6.037e+04 1.404e+04 4.299 6.07e-05 ***
M7 5.658e+04 1.245e+04 4.545 2.54e-05 ***
M8 3.692e+04 1.381e+04 2.674 0.009548 **
M9 1.135e+05 1.455e+04 7.796 8.10e-11 ***
M10 8.686e+04 1.396e+04 6.223 4.42e-08 ***
M11 7.303e+04 1.413e+04 5.170 2.58e-06 ***
t 1.263e+02 1.114e+02 1.133 0.261431
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 20440 on 63 degrees of freedom
Multiple R-squared: 0.8106, Adjusted R-squared: 0.7625
F-statistic: 16.85 on 16 and 63 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.1734070770 0.3468141540 0.8265929
[2,] 0.0905598261 0.1811196522 0.9094402
[3,] 0.0422410208 0.0844820417 0.9577590
[4,] 0.0176546392 0.0353092783 0.9823454
[5,] 0.0069317753 0.0138635505 0.9930682
[6,] 0.1531938519 0.3063877039 0.8468061
[7,] 0.1291532779 0.2583065558 0.8708467
[8,] 0.0794178600 0.1588357200 0.9205821
[9,] 0.1347318639 0.2694637278 0.8652681
[10,] 0.0891623579 0.1783247159 0.9108376
[11,] 0.0586686663 0.1173373326 0.9413313
[12,] 0.0423934297 0.0847868594 0.9576066
[13,] 0.0248223612 0.0496447223 0.9751776
[14,] 0.0161014635 0.0322029270 0.9838985
[15,] 0.0093334201 0.0186668402 0.9906666
[16,] 0.0070668792 0.0141337585 0.9929331
[17,] 0.0041210430 0.0082420859 0.9958790
[18,] 0.0075554985 0.0151109969 0.9924445
[19,] 0.0040018819 0.0080037638 0.9959981
[20,] 0.0032477548 0.0064955096 0.9967522
[21,] 0.0024196047 0.0048392094 0.9975804
[22,] 0.0016667359 0.0033334717 0.9983333
[23,] 0.0012378936 0.0024757873 0.9987621
[24,] 0.0006782197 0.0013564395 0.9993218
[25,] 0.0005060979 0.0010121958 0.9994939
[26,] 0.0003890860 0.0007781720 0.9996109
[27,] 0.0001989038 0.0003978075 0.9998011
[28,] 0.0104754695 0.0209509391 0.9895245
[29,] 0.0847344297 0.1694688594 0.9152656
[30,] 0.0684975989 0.1369951977 0.9315024
[31,] 0.0535883755 0.1071767510 0.9464116
[32,] 0.0521103493 0.1042206987 0.9478897
[33,] 0.0318213928 0.0636427855 0.9681786
[34,] 0.0347253906 0.0694507812 0.9652746
[35,] 0.0216942254 0.0433884508 0.9783058
[36,] 0.0123140847 0.0246281695 0.9876859
[37,] 0.0764534804 0.1529069608 0.9235465
[38,] 0.0530612027 0.1061224055 0.9469388
[39,] 0.0446639805 0.0893279609 0.9553360
[40,] 0.1732372493 0.3464744985 0.8267628
[41,] 0.1007720001 0.2015440003 0.8992280
> postscript(file="/var/www/html/rcomp/tmp/1coqn1293184469.ps",horizontal=F,onefile=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/2coqn1293184469.ps",horizontal=F,onefile=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/3nfp81293184469.ps",horizontal=F,onefile=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/4nfp81293184469.ps",horizontal=F,onefile=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/5nfp81293184469.ps",horizontal=F,onefile=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 = 80
Frequency = 1
1 2 3 4 5 6
-20321.2658 25638.7066 7378.3782 9664.7981 3935.5118 11411.3463
7 8 9 10 11 12
8769.7473 -4906.2043 -4185.6134 -5282.7713 -4863.2008 -2979.0671
13 14 15 16 17 18
5869.7545 -5636.7958 -4814.2734 17647.1587 -6019.2673 6366.5658
19 20 21 22 23 24
-669.7719 5508.0002 7638.6992 2312.1762 -10886.5189 2479.4678
25 26 27 28 29 30
51497.3959 -23129.4904 -5934.5680 -10731.9980 5878.5569 -2263.7826
31 32 33 34 35 36
-570.7232 -1591.1548 11633.1524 1860.8259 -25939.5234 -7428.8544
37 38 39 40 41 42
-12376.8584 -2494.6862 16896.7261 -12477.3632 9019.6671 -8956.3066
43 44 45 46 47 48
-7501.5383 7607.5229 -7682.6584 -15478.6868 -54659.4721 35656.0923
49 50 51 52 53 54
-27394.2883 10030.2424 -106.2312 -6227.7314 6756.0862 -19852.3571
55 56 57 58 59 60
-10096.1047 -3394.5993 -23677.8893 -18320.3436 6534.5964 -5558.5205
61 62 63 64 65 66
-9728.1966 -9014.6314 7125.0851 7018.2558 -5264.0110 -6987.0537
67 68 69 70 71 72
12786.5942 7535.2981 16274.3094 34908.7997 89814.1189 -22169.1181
73 74 75 76 77 78
12453.4586 4606.6547 -20545.1168 -4893.1200 -14306.5437 20281.5879
79 80
-2718.2035 -10758.8628
> postscript(file="/var/www/html/rcomp/tmp/6xopt1293184469.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> dum <- cbind(lag(myerror,k=1),myerror)
> dum
Time Series:
Start = 0
End = 80
Frequency = 1
lag(myerror, k = 1) myerror
0 -20321.2658 NA
1 25638.7066 -20321.2658
2 7378.3782 25638.7066
3 9664.7981 7378.3782
4 3935.5118 9664.7981
5 11411.3463 3935.5118
6 8769.7473 11411.3463
7 -4906.2043 8769.7473
8 -4185.6134 -4906.2043
9 -5282.7713 -4185.6134
10 -4863.2008 -5282.7713
11 -2979.0671 -4863.2008
12 5869.7545 -2979.0671
13 -5636.7958 5869.7545
14 -4814.2734 -5636.7958
15 17647.1587 -4814.2734
16 -6019.2673 17647.1587
17 6366.5658 -6019.2673
18 -669.7719 6366.5658
19 5508.0002 -669.7719
20 7638.6992 5508.0002
21 2312.1762 7638.6992
22 -10886.5189 2312.1762
23 2479.4678 -10886.5189
24 51497.3959 2479.4678
25 -23129.4904 51497.3959
26 -5934.5680 -23129.4904
27 -10731.9980 -5934.5680
28 5878.5569 -10731.9980
29 -2263.7826 5878.5569
30 -570.7232 -2263.7826
31 -1591.1548 -570.7232
32 11633.1524 -1591.1548
33 1860.8259 11633.1524
34 -25939.5234 1860.8259
35 -7428.8544 -25939.5234
36 -12376.8584 -7428.8544
37 -2494.6862 -12376.8584
38 16896.7261 -2494.6862
39 -12477.3632 16896.7261
40 9019.6671 -12477.3632
41 -8956.3066 9019.6671
42 -7501.5383 -8956.3066
43 7607.5229 -7501.5383
44 -7682.6584 7607.5229
45 -15478.6868 -7682.6584
46 -54659.4721 -15478.6868
47 35656.0923 -54659.4721
48 -27394.2883 35656.0923
49 10030.2424 -27394.2883
50 -106.2312 10030.2424
51 -6227.7314 -106.2312
52 6756.0862 -6227.7314
53 -19852.3571 6756.0862
54 -10096.1047 -19852.3571
55 -3394.5993 -10096.1047
56 -23677.8893 -3394.5993
57 -18320.3436 -23677.8893
58 6534.5964 -18320.3436
59 -5558.5205 6534.5964
60 -9728.1966 -5558.5205
61 -9014.6314 -9728.1966
62 7125.0851 -9014.6314
63 7018.2558 7125.0851
64 -5264.0110 7018.2558
65 -6987.0537 -5264.0110
66 12786.5942 -6987.0537
67 7535.2981 12786.5942
68 16274.3094 7535.2981
69 34908.7997 16274.3094
70 89814.1189 34908.7997
71 -22169.1181 89814.1189
72 12453.4586 -22169.1181
73 4606.6547 12453.4586
74 -20545.1168 4606.6547
75 -4893.1200 -20545.1168
76 -14306.5437 -4893.1200
77 20281.5879 -14306.5437
78 -2718.2035 20281.5879
79 -10758.8628 -2718.2035
80 NA -10758.8628
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 25638.7066 -20321.2658
[2,] 7378.3782 25638.7066
[3,] 9664.7981 7378.3782
[4,] 3935.5118 9664.7981
[5,] 11411.3463 3935.5118
[6,] 8769.7473 11411.3463
[7,] -4906.2043 8769.7473
[8,] -4185.6134 -4906.2043
[9,] -5282.7713 -4185.6134
[10,] -4863.2008 -5282.7713
[11,] -2979.0671 -4863.2008
[12,] 5869.7545 -2979.0671
[13,] -5636.7958 5869.7545
[14,] -4814.2734 -5636.7958
[15,] 17647.1587 -4814.2734
[16,] -6019.2673 17647.1587
[17,] 6366.5658 -6019.2673
[18,] -669.7719 6366.5658
[19,] 5508.0002 -669.7719
[20,] 7638.6992 5508.0002
[21,] 2312.1762 7638.6992
[22,] -10886.5189 2312.1762
[23,] 2479.4678 -10886.5189
[24,] 51497.3959 2479.4678
[25,] -23129.4904 51497.3959
[26,] -5934.5680 -23129.4904
[27,] -10731.9980 -5934.5680
[28,] 5878.5569 -10731.9980
[29,] -2263.7826 5878.5569
[30,] -570.7232 -2263.7826
[31,] -1591.1548 -570.7232
[32,] 11633.1524 -1591.1548
[33,] 1860.8259 11633.1524
[34,] -25939.5234 1860.8259
[35,] -7428.8544 -25939.5234
[36,] -12376.8584 -7428.8544
[37,] -2494.6862 -12376.8584
[38,] 16896.7261 -2494.6862
[39,] -12477.3632 16896.7261
[40,] 9019.6671 -12477.3632
[41,] -8956.3066 9019.6671
[42,] -7501.5383 -8956.3066
[43,] 7607.5229 -7501.5383
[44,] -7682.6584 7607.5229
[45,] -15478.6868 -7682.6584
[46,] -54659.4721 -15478.6868
[47,] 35656.0923 -54659.4721
[48,] -27394.2883 35656.0923
[49,] 10030.2424 -27394.2883
[50,] -106.2312 10030.2424
[51,] -6227.7314 -106.2312
[52,] 6756.0862 -6227.7314
[53,] -19852.3571 6756.0862
[54,] -10096.1047 -19852.3571
[55,] -3394.5993 -10096.1047
[56,] -23677.8893 -3394.5993
[57,] -18320.3436 -23677.8893
[58,] 6534.5964 -18320.3436
[59,] -5558.5205 6534.5964
[60,] -9728.1966 -5558.5205
[61,] -9014.6314 -9728.1966
[62,] 7125.0851 -9014.6314
[63,] 7018.2558 7125.0851
[64,] -5264.0110 7018.2558
[65,] -6987.0537 -5264.0110
[66,] 12786.5942 -6987.0537
[67,] 7535.2981 12786.5942
[68,] 16274.3094 7535.2981
[69,] 34908.7997 16274.3094
[70,] 89814.1189 34908.7997
[71,] -22169.1181 89814.1189
[72,] 12453.4586 -22169.1181
[73,] 4606.6547 12453.4586
[74,] -20545.1168 4606.6547
[75,] -4893.1200 -20545.1168
[76,] -14306.5437 -4893.1200
[77,] 20281.5879 -14306.5437
[78,] -2718.2035 20281.5879
[79,] -10758.8628 -2718.2035
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 25638.7066 -20321.2658
2 7378.3782 25638.7066
3 9664.7981 7378.3782
4 3935.5118 9664.7981
5 11411.3463 3935.5118
6 8769.7473 11411.3463
7 -4906.2043 8769.7473
8 -4185.6134 -4906.2043
9 -5282.7713 -4185.6134
10 -4863.2008 -5282.7713
11 -2979.0671 -4863.2008
12 5869.7545 -2979.0671
13 -5636.7958 5869.7545
14 -4814.2734 -5636.7958
15 17647.1587 -4814.2734
16 -6019.2673 17647.1587
17 6366.5658 -6019.2673
18 -669.7719 6366.5658
19 5508.0002 -669.7719
20 7638.6992 5508.0002
21 2312.1762 7638.6992
22 -10886.5189 2312.1762
23 2479.4678 -10886.5189
24 51497.3959 2479.4678
25 -23129.4904 51497.3959
26 -5934.5680 -23129.4904
27 -10731.9980 -5934.5680
28 5878.5569 -10731.9980
29 -2263.7826 5878.5569
30 -570.7232 -2263.7826
31 -1591.1548 -570.7232
32 11633.1524 -1591.1548
33 1860.8259 11633.1524
34 -25939.5234 1860.8259
35 -7428.8544 -25939.5234
36 -12376.8584 -7428.8544
37 -2494.6862 -12376.8584
38 16896.7261 -2494.6862
39 -12477.3632 16896.7261
40 9019.6671 -12477.3632
41 -8956.3066 9019.6671
42 -7501.5383 -8956.3066
43 7607.5229 -7501.5383
44 -7682.6584 7607.5229
45 -15478.6868 -7682.6584
46 -54659.4721 -15478.6868
47 35656.0923 -54659.4721
48 -27394.2883 35656.0923
49 10030.2424 -27394.2883
50 -106.2312 10030.2424
51 -6227.7314 -106.2312
52 6756.0862 -6227.7314
53 -19852.3571 6756.0862
54 -10096.1047 -19852.3571
55 -3394.5993 -10096.1047
56 -23677.8893 -3394.5993
57 -18320.3436 -23677.8893
58 6534.5964 -18320.3436
59 -5558.5205 6534.5964
60 -9728.1966 -5558.5205
61 -9014.6314 -9728.1966
62 7125.0851 -9014.6314
63 7018.2558 7125.0851
64 -5264.0110 7018.2558
65 -6987.0537 -5264.0110
66 12786.5942 -6987.0537
67 7535.2981 12786.5942
68 16274.3094 7535.2981
69 34908.7997 16274.3094
70 89814.1189 34908.7997
71 -22169.1181 89814.1189
72 12453.4586 -22169.1181
73 4606.6547 12453.4586
74 -20545.1168 4606.6547
75 -4893.1200 -20545.1168
76 -14306.5437 -4893.1200
77 20281.5879 -14306.5437
78 -2718.2035 20281.5879
79 -10758.8628 -2718.2035
> 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/78foe1293184469.ps",horizontal=F,onefile=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/88foe1293184469.ps",horizontal=F,onefile=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/98foe1293184469.ps",horizontal=F,onefile=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/10jp5h1293184469.ps",horizontal=F,onefile=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/11mpmn1293184469.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/12882s1293184469.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/134iij1293184469.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/14w9zm1293184469.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/1509ga1293184469.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/16e1ej1293184469.tab")
+ }
>
> try(system("convert tmp/1coqn1293184469.ps tmp/1coqn1293184469.png",intern=TRUE))
character(0)
> try(system("convert tmp/2coqn1293184469.ps tmp/2coqn1293184469.png",intern=TRUE))
character(0)
> try(system("convert tmp/3nfp81293184469.ps tmp/3nfp81293184469.png",intern=TRUE))
character(0)
> try(system("convert tmp/4nfp81293184469.ps tmp/4nfp81293184469.png",intern=TRUE))
character(0)
> try(system("convert tmp/5nfp81293184469.ps tmp/5nfp81293184469.png",intern=TRUE))
character(0)
> try(system("convert tmp/6xopt1293184469.ps tmp/6xopt1293184469.png",intern=TRUE))
character(0)
> try(system("convert tmp/78foe1293184469.ps tmp/78foe1293184469.png",intern=TRUE))
character(0)
> try(system("convert tmp/88foe1293184469.ps tmp/88foe1293184469.png",intern=TRUE))
character(0)
> try(system("convert tmp/98foe1293184469.ps tmp/98foe1293184469.png",intern=TRUE))
character(0)
> try(system("convert tmp/10jp5h1293184469.ps tmp/10jp5h1293184469.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
2.752 1.678 8.506