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(264768
+ ,236089
+ ,268586
+ ,268361
+ ,272433
+ ,274412
+ ,269974
+ ,236997
+ ,264768
+ ,268586
+ ,268361
+ ,272433
+ ,304744
+ ,264579
+ ,269974
+ ,264768
+ ,268586
+ ,268361
+ ,309365
+ ,270349
+ ,304744
+ ,269974
+ ,264768
+ ,268586
+ ,308347
+ ,269645
+ ,309365
+ ,304744
+ ,269974
+ ,264768
+ ,298427
+ ,267037
+ ,308347
+ ,309365
+ ,304744
+ ,269974
+ ,289231
+ ,258113
+ ,298427
+ ,308347
+ ,309365
+ ,304744
+ ,291975
+ ,262813
+ ,289231
+ ,298427
+ ,308347
+ ,309365
+ ,294912
+ ,267413
+ ,291975
+ ,289231
+ ,298427
+ ,308347
+ ,293488
+ ,267366
+ ,294912
+ ,291975
+ ,289231
+ ,298427
+ ,290555
+ ,264777
+ ,293488
+ ,294912
+ ,291975
+ ,289231
+ ,284736
+ ,258863
+ ,290555
+ ,293488
+ ,294912
+ ,291975
+ ,281818
+ ,254844
+ ,284736
+ ,290555
+ ,293488
+ ,294912
+ ,287854
+ ,254868
+ ,281818
+ ,284736
+ ,290555
+ ,293488
+ ,316263
+ ,277267
+ ,287854
+ ,281818
+ ,284736
+ ,290555
+ ,325412
+ ,285351
+ ,316263
+ ,287854
+ ,281818
+ ,284736
+ ,326011
+ ,286602
+ ,325412
+ ,316263
+ ,287854
+ ,281818
+ ,328282
+ ,283042
+ ,326011
+ ,325412
+ ,316263
+ ,287854
+ ,317480
+ ,276687
+ ,328282
+ ,326011
+ ,325412
+ ,316263
+ ,317539
+ ,277915
+ ,317480
+ ,328282
+ ,326011
+ ,325412
+ ,313737
+ ,277128
+ ,317539
+ ,317480
+ ,328282
+ ,326011
+ ,312276
+ ,277103
+ ,313737
+ ,317539
+ ,317480
+ ,328282
+ ,309391
+ ,275037
+ ,312276
+ ,313737
+ ,317539
+ ,317480
+ ,302950
+ ,270150
+ ,309391
+ ,312276
+ ,313737
+ ,317539
+ ,300316
+ ,267140
+ ,302950
+ ,309391
+ ,312276
+ ,313737
+ ,304035
+ ,264993
+ ,300316
+ ,302950
+ ,309391
+ ,312276
+ ,333476
+ ,287259
+ ,304035
+ ,300316
+ ,302950
+ ,309391
+ ,337698
+ ,291186
+ ,333476
+ ,304035
+ ,300316
+ ,302950
+ ,335932
+ ,292300
+ ,337698
+ ,333476
+ ,304035
+ ,300316
+ ,323931
+ ,288186
+ ,335932
+ ,337698
+ ,333476
+ ,304035
+ ,313927
+ ,281477
+ ,323931
+ ,335932
+ ,337698
+ ,333476
+ ,314485
+ ,282656
+ ,313927
+ ,323931
+ ,335932
+ ,337698
+ ,313218
+ ,280190
+ ,314485
+ ,313927
+ ,323931
+ ,335932
+ ,309664
+ ,280408
+ ,313218
+ ,314485
+ ,313927
+ ,323931
+ ,302963
+ ,276836
+ ,309664
+ ,313218
+ ,314485
+ ,313927
+ ,298989
+ ,275216
+ ,302963
+ ,309664
+ ,313218
+ ,314485
+ ,298423
+ ,274352
+ ,298989
+ ,302963
+ ,309664
+ ,313218
+ ,301631
+ ,271311
+ ,298423
+ ,298989
+ ,302963
+ ,309664
+ ,329765
+ ,289802
+ ,301631
+ ,298423
+ ,298989
+ ,302963
+ ,335083
+ ,290726
+ ,329765
+ ,301631
+ ,298423
+ ,298989
+ ,327616
+ ,292300
+ ,335083
+ ,329765
+ ,301631
+ ,298423
+ ,309119
+ ,278506
+ ,327616
+ ,335083
+ ,329765
+ ,301631
+ ,295916
+ ,269826
+ ,309119
+ ,327616
+ ,335083
+ ,329765
+ ,291413
+ ,265861
+ ,295916
+ ,309119
+ ,327616
+ ,335083
+ ,291542
+ ,269034
+ ,291413
+ ,295916
+ ,309119
+ ,327616
+ ,284678
+ ,264176
+ ,291542
+ ,291413
+ ,295916
+ ,309119
+ ,276475
+ ,255198
+ ,284678
+ ,291542
+ ,291413
+ ,295916
+ ,272566
+ ,253353
+ ,276475
+ ,284678
+ ,291542
+ ,291413
+ ,264981
+ ,246057
+ ,272566
+ ,276475
+ ,284678
+ ,291542
+ ,263290
+ ,235372
+ ,264981
+ ,272566
+ ,276475
+ ,284678
+ ,296806
+ ,258556
+ ,263290
+ ,264981
+ ,272566
+ ,276475
+ ,303598
+ ,260993
+ ,296806
+ ,263290
+ ,264981
+ ,272566
+ ,286994
+ ,254663
+ ,303598
+ ,296806
+ ,263290
+ ,264981
+ ,276427
+ ,250643
+ ,286994
+ ,303598
+ ,296806
+ ,263290
+ ,266424
+ ,243422
+ ,276427
+ ,286994
+ ,303598
+ ,296806
+ ,267153
+ ,247105
+ ,266424
+ ,276427
+ ,286994
+ ,303598
+ ,268381
+ ,248541
+ ,267153
+ ,266424
+ ,276427
+ ,286994
+ ,262522
+ ,245039
+ ,268381
+ ,267153
+ ,266424
+ ,276427
+ ,255542
+ ,237080
+ ,262522
+ ,268381
+ ,267153
+ ,266424
+ ,253158
+ ,237085
+ ,255542
+ ,262522
+ ,268381
+ ,267153
+ ,243803
+ ,225554
+ ,253158
+ ,255542
+ ,262522
+ ,268381
+ ,250741
+ ,226839
+ ,243803
+ ,253158
+ ,255542
+ ,262522
+ ,280445
+ ,247934
+ ,250741
+ ,243803
+ ,253158
+ ,255542
+ ,285257
+ ,248333
+ ,280445
+ ,250741
+ ,243803
+ ,253158
+ ,270976
+ ,246969
+ ,285257
+ ,280445
+ ,250741
+ ,243803
+ ,261076
+ ,245098
+ ,270976
+ ,285257
+ ,280445
+ ,250741
+ ,255603
+ ,246263
+ ,261076
+ ,270976
+ ,285257
+ ,280445
+ ,260376
+ ,255765
+ ,255603
+ ,261076
+ ,270976
+ ,285257
+ ,263903
+ ,264319
+ ,260376
+ ,255603
+ ,261076
+ ,270976
+ ,264291
+ ,268347
+ ,263903
+ ,260376
+ ,255603
+ ,261076
+ ,263276
+ ,273046
+ ,264291
+ ,263903
+ ,260376
+ ,255603
+ ,262572
+ ,273963
+ ,263276
+ ,264291
+ ,263903
+ ,260376
+ ,256167
+ ,267430
+ ,262572
+ ,263276
+ ,264291
+ ,263903
+ ,264221
+ ,271993
+ ,256167
+ ,262572
+ ,263276
+ ,264291
+ ,293860
+ ,292710
+ ,264221
+ ,256167
+ ,262572
+ ,263276)
+ ,dim=c(6
+ ,75)
+ ,dimnames=list(c('y'
+ ,'x'
+ ,'y1'
+ ,'y2'
+ ,'y3'
+ ,'y4')
+ ,1:75))
> y <- array(NA,dim=c(6,75),dimnames=list(c('y','x','y1','y2','y3','y4'),1:75))
> 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 = '2'
> #'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
x y y1 y2 y3 y4 M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 236089 264768 268586 268361 272433 274412 1 0 0 0 0 0 0 0 0 0 0
2 236997 269974 264768 268586 268361 272433 0 1 0 0 0 0 0 0 0 0 0
3 264579 304744 269974 264768 268586 268361 0 0 1 0 0 0 0 0 0 0 0
4 270349 309365 304744 269974 264768 268586 0 0 0 1 0 0 0 0 0 0 0
5 269645 308347 309365 304744 269974 264768 0 0 0 0 1 0 0 0 0 0 0
6 267037 298427 308347 309365 304744 269974 0 0 0 0 0 1 0 0 0 0 0
7 258113 289231 298427 308347 309365 304744 0 0 0 0 0 0 1 0 0 0 0
8 262813 291975 289231 298427 308347 309365 0 0 0 0 0 0 0 1 0 0 0
9 267413 294912 291975 289231 298427 308347 0 0 0 0 0 0 0 0 1 0 0
10 267366 293488 294912 291975 289231 298427 0 0 0 0 0 0 0 0 0 1 0
11 264777 290555 293488 294912 291975 289231 0 0 0 0 0 0 0 0 0 0 1
12 258863 284736 290555 293488 294912 291975 0 0 0 0 0 0 0 0 0 0 0
13 254844 281818 284736 290555 293488 294912 1 0 0 0 0 0 0 0 0 0 0
14 254868 287854 281818 284736 290555 293488 0 1 0 0 0 0 0 0 0 0 0
15 277267 316263 287854 281818 284736 290555 0 0 1 0 0 0 0 0 0 0 0
16 285351 325412 316263 287854 281818 284736 0 0 0 1 0 0 0 0 0 0 0
17 286602 326011 325412 316263 287854 281818 0 0 0 0 1 0 0 0 0 0 0
18 283042 328282 326011 325412 316263 287854 0 0 0 0 0 1 0 0 0 0 0
19 276687 317480 328282 326011 325412 316263 0 0 0 0 0 0 1 0 0 0 0
20 277915 317539 317480 328282 326011 325412 0 0 0 0 0 0 0 1 0 0 0
21 277128 313737 317539 317480 328282 326011 0 0 0 0 0 0 0 0 1 0 0
22 277103 312276 313737 317539 317480 328282 0 0 0 0 0 0 0 0 0 1 0
23 275037 309391 312276 313737 317539 317480 0 0 0 0 0 0 0 0 0 0 1
24 270150 302950 309391 312276 313737 317539 0 0 0 0 0 0 0 0 0 0 0
25 267140 300316 302950 309391 312276 313737 1 0 0 0 0 0 0 0 0 0 0
26 264993 304035 300316 302950 309391 312276 0 1 0 0 0 0 0 0 0 0 0
27 287259 333476 304035 300316 302950 309391 0 0 1 0 0 0 0 0 0 0 0
28 291186 337698 333476 304035 300316 302950 0 0 0 1 0 0 0 0 0 0 0
29 292300 335932 337698 333476 304035 300316 0 0 0 0 1 0 0 0 0 0 0
30 288186 323931 335932 337698 333476 304035 0 0 0 0 0 1 0 0 0 0 0
31 281477 313927 323931 335932 337698 333476 0 0 0 0 0 0 1 0 0 0 0
32 282656 314485 313927 323931 335932 337698 0 0 0 0 0 0 0 1 0 0 0
33 280190 313218 314485 313927 323931 335932 0 0 0 0 0 0 0 0 1 0 0
34 280408 309664 313218 314485 313927 323931 0 0 0 0 0 0 0 0 0 1 0
35 276836 302963 309664 313218 314485 313927 0 0 0 0 0 0 0 0 0 0 1
36 275216 298989 302963 309664 313218 314485 0 0 0 0 0 0 0 0 0 0 0
37 274352 298423 298989 302963 309664 313218 1 0 0 0 0 0 0 0 0 0 0
38 271311 301631 298423 298989 302963 309664 0 1 0 0 0 0 0 0 0 0 0
39 289802 329765 301631 298423 298989 302963 0 0 1 0 0 0 0 0 0 0 0
40 290726 335083 329765 301631 298423 298989 0 0 0 1 0 0 0 0 0 0 0
41 292300 327616 335083 329765 301631 298423 0 0 0 0 1 0 0 0 0 0 0
42 278506 309119 327616 335083 329765 301631 0 0 0 0 0 1 0 0 0 0 0
43 269826 295916 309119 327616 335083 329765 0 0 0 0 0 0 1 0 0 0 0
44 265861 291413 295916 309119 327616 335083 0 0 0 0 0 0 0 1 0 0 0
45 269034 291542 291413 295916 309119 327616 0 0 0 0 0 0 0 0 1 0 0
46 264176 284678 291542 291413 295916 309119 0 0 0 0 0 0 0 0 0 1 0
47 255198 276475 284678 291542 291413 295916 0 0 0 0 0 0 0 0 0 0 1
48 253353 272566 276475 284678 291542 291413 0 0 0 0 0 0 0 0 0 0 0
49 246057 264981 272566 276475 284678 291542 1 0 0 0 0 0 0 0 0 0 0
50 235372 263290 264981 272566 276475 284678 0 1 0 0 0 0 0 0 0 0 0
51 258556 296806 263290 264981 272566 276475 0 0 1 0 0 0 0 0 0 0 0
52 260993 303598 296806 263290 264981 272566 0 0 0 1 0 0 0 0 0 0 0
53 254663 286994 303598 296806 263290 264981 0 0 0 0 1 0 0 0 0 0 0
54 250643 276427 286994 303598 296806 263290 0 0 0 0 0 1 0 0 0 0 0
55 243422 266424 276427 286994 303598 296806 0 0 0 0 0 0 1 0 0 0 0
56 247105 267153 266424 276427 286994 303598 0 0 0 0 0 0 0 1 0 0 0
57 248541 268381 267153 266424 276427 286994 0 0 0 0 0 0 0 0 1 0 0
58 245039 262522 268381 267153 266424 276427 0 0 0 0 0 0 0 0 0 1 0
59 237080 255542 262522 268381 267153 266424 0 0 0 0 0 0 0 0 0 0 1
60 237085 253158 255542 262522 268381 267153 0 0 0 0 0 0 0 0 0 0 0
61 225554 243803 253158 255542 262522 268381 1 0 0 0 0 0 0 0 0 0 0
62 226839 250741 243803 253158 255542 262522 0 1 0 0 0 0 0 0 0 0 0
63 247934 280445 250741 243803 253158 255542 0 0 1 0 0 0 0 0 0 0 0
64 248333 285257 280445 250741 243803 253158 0 0 0 1 0 0 0 0 0 0 0
65 246969 270976 285257 280445 250741 243803 0 0 0 0 1 0 0 0 0 0 0
66 245098 261076 270976 285257 280445 250741 0 0 0 0 0 1 0 0 0 0 0
67 246263 255603 261076 270976 285257 280445 0 0 0 0 0 0 1 0 0 0 0
68 255765 260376 255603 261076 270976 285257 0 0 0 0 0 0 0 1 0 0 0
69 264319 263903 260376 255603 261076 270976 0 0 0 0 0 0 0 0 1 0 0
70 268347 264291 263903 260376 255603 261076 0 0 0 0 0 0 0 0 0 1 0
71 273046 263276 264291 263903 260376 255603 0 0 0 0 0 0 0 0 0 0 1
72 273963 262572 263276 264291 263903 260376 0 0 0 0 0 0 0 0 0 0 0
73 267430 256167 262572 263276 264291 263903 1 0 0 0 0 0 0 0 0 0 0
74 271993 264221 256167 262572 263276 264291 0 1 0 0 0 0 0 0 0 0 0
75 292710 293860 264221 256167 262572 263276 0 0 1 0 0 0 0 0 0 0 0
t
1 1
2 2
3 3
4 4
5 5
6 6
7 7
8 8
9 9
10 10
11 11
12 12
13 13
14 14
15 15
16 16
17 17
18 18
19 19
20 20
21 21
22 22
23 23
24 24
25 25
26 26
27 27
28 28
29 29
30 30
31 31
32 32
33 33
34 34
35 35
36 36
37 37
38 38
39 39
40 40
41 41
42 42
43 43
44 44
45 45
46 46
47 47
48 48
49 49
50 50
51 51
52 52
53 53
54 54
55 55
56 56
57 57
58 58
59 59
60 60
61 61
62 62
63 63
64 64
65 65
66 66
67 67
68 68
69 69
70 70
71 71
72 72
73 73
74 74
75 75
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) y y1 y2 y3 y4
2.941e+04 1.281e+00 2.137e-02 -1.100e-01 1.891e-01 -5.835e-01
M1 M2 M3 M4 M5 M6
6.450e+02 -7.863e+03 -2.777e+04 -3.605e+04 -2.860e+04 -2.417e+04
M7 M8 M9 M10 M11 t
-1.856e+03 3.395e+03 1.693e+03 8.740e+02 -2.482e+03 3.477e+02
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-12215.5 -2825.5 475.4 2773.1 15156.3
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.941e+04 1.568e+04 1.875 0.065860 .
y 1.281e+00 2.402e-01 5.332 1.73e-06 ***
y1 2.137e-02 3.549e-01 0.060 0.952206
y2 -1.100e-01 3.574e-01 -0.308 0.759483
y3 1.891e-01 3.606e-01 0.524 0.601970
y4 -5.835e-01 2.364e-01 -2.468 0.016614 *
M1 6.450e+02 3.895e+03 0.166 0.869045
M2 -7.863e+03 4.443e+03 -1.770 0.082125 .
M3 -2.777e+04 9.119e+03 -3.045 0.003515 **
M4 -3.605e+04 9.352e+03 -3.855 0.000296 ***
M5 -2.860e+04 8.985e+03 -3.183 0.002362 **
M6 -2.417e+04 8.377e+03 -2.885 0.005517 **
M7 -1.856e+03 4.558e+03 -0.407 0.685340
M8 3.395e+03 4.771e+03 0.711 0.479682
M9 1.693e+03 5.147e+03 0.329 0.743500
M10 8.740e+02 4.915e+03 0.178 0.859490
M11 -2.482e+03 3.989e+03 -0.622 0.536261
t 3.477e+02 4.925e+01 7.060 2.53e-09 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 6717 on 57 degrees of freedom
Multiple R-squared: 0.8762, Adjusted R-squared: 0.8393
F-statistic: 23.73 on 17 and 57 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.34992333 0.699846668 0.650076666
[2,] 0.23606819 0.472136378 0.763931811
[3,] 0.14606601 0.292132016 0.853933992
[4,] 0.08884115 0.177682297 0.911158851
[5,] 0.04683939 0.093678779 0.953160610
[6,] 0.04034923 0.080698454 0.959650773
[7,] 0.07693181 0.153863629 0.923068185
[8,] 0.62229992 0.755400164 0.377700082
[9,] 0.61461652 0.770766956 0.385383478
[10,] 0.55736570 0.885268597 0.442634298
[11,] 0.48141524 0.962830479 0.518584760
[12,] 0.42817326 0.856346514 0.571826743
[13,] 0.46805586 0.936111727 0.531944137
[14,] 0.42275889 0.845517771 0.577241115
[15,] 0.34289957 0.685799134 0.657100433
[16,] 0.28288548 0.565770959 0.717114521
[17,] 0.28016071 0.560321415 0.719839293
[18,] 0.23093551 0.461871020 0.769064490
[19,] 0.18163188 0.363263758 0.818368121
[20,] 0.28882212 0.577644250 0.711177875
[21,] 0.50770087 0.984598265 0.492299133
[22,] 0.88340092 0.233198160 0.116599080
[23,] 0.91996691 0.160066183 0.080033091
[24,] 0.94149742 0.117005168 0.058502584
[25,] 0.92346569 0.153068611 0.076534305
[26,] 0.90359104 0.192817910 0.096408955
[27,] 0.89673572 0.206528568 0.103264284
[28,] 0.86025671 0.279486576 0.139743288
[29,] 0.83761754 0.324764930 0.162382465
[30,] 0.99753910 0.004921797 0.002460898
[31,] 0.99732242 0.005355169 0.002677584
[32,] 0.99388037 0.012239250 0.006119625
[33,] 0.97951327 0.040973456 0.020486728
[34,] 0.94204306 0.115913888 0.057956944
> postscript(file="/var/www/html/rcomp/tmp/1dson1258904733.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/2d3sm1258904733.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/3ymb91258904733.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/4wc8n1258904733.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/59q9s1258904733.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 = 75
Frequency = 1
1 2 3 4 5 6
-1015.7160 1107.1534 774.8153 9245.3258 2552.8818 4861.9943
7 8 9 10 11 12
4567.9963 2150.1252 4555.7867 2993.3220 1637.3861 1297.0945
13 14 15 16 17 18
1806.9008 1407.5119 5926.8749 7443.3713 207.2118 -11896.6661
19 20 21 22 23 24
-12215.5416 -10956.1878 -6789.0778 -1016.6045 -3080.4166 -1894.3420
25 26 27 28 29 30
-4645.6402 -4354.0344 -1062.8592 1912.6503 -1609.2611 1971.1196
31 32 33 34 35 36
1855.2232 -1587.1649 -948.9054 -730.6276 1280.9572 2237.9376
37 38 39 40 41 42
387.0364 166.6083 -1098.2838 -1509.2634 3865.5799 6276.0097
43 44 45 46 47 48
6827.7963 5794.2706 7942.4728 3552.0242 1395.9397 -1504.1158
49 50 51 52 53 54
475.4229 -2605.4953 -7625.2404 -7699.0269 -1136.0075 -2625.6339
55 56 57 58 59 60
-3025.3468 280.7475 -7306.2939 -7054.2695 -8780.5795 -8854.3142
61 62 63 64 65 66
-8290.0955 -9890.6219 -12071.5997 -9393.0571 -3880.4048 1413.1764
67 68 69 70 71 72
1989.8726 4318.2094 2546.0176 2256.1555 7546.7131 8717.7399
73 74 75
11282.0916 14168.8779 15156.2927
> postscript(file="/var/www/html/rcomp/tmp/68hnx1258904733.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 = 75
Frequency = 1
lag(myerror, k = 1) myerror
0 -1015.7160 NA
1 1107.1534 -1015.7160
2 774.8153 1107.1534
3 9245.3258 774.8153
4 2552.8818 9245.3258
5 4861.9943 2552.8818
6 4567.9963 4861.9943
7 2150.1252 4567.9963
8 4555.7867 2150.1252
9 2993.3220 4555.7867
10 1637.3861 2993.3220
11 1297.0945 1637.3861
12 1806.9008 1297.0945
13 1407.5119 1806.9008
14 5926.8749 1407.5119
15 7443.3713 5926.8749
16 207.2118 7443.3713
17 -11896.6661 207.2118
18 -12215.5416 -11896.6661
19 -10956.1878 -12215.5416
20 -6789.0778 -10956.1878
21 -1016.6045 -6789.0778
22 -3080.4166 -1016.6045
23 -1894.3420 -3080.4166
24 -4645.6402 -1894.3420
25 -4354.0344 -4645.6402
26 -1062.8592 -4354.0344
27 1912.6503 -1062.8592
28 -1609.2611 1912.6503
29 1971.1196 -1609.2611
30 1855.2232 1971.1196
31 -1587.1649 1855.2232
32 -948.9054 -1587.1649
33 -730.6276 -948.9054
34 1280.9572 -730.6276
35 2237.9376 1280.9572
36 387.0364 2237.9376
37 166.6083 387.0364
38 -1098.2838 166.6083
39 -1509.2634 -1098.2838
40 3865.5799 -1509.2634
41 6276.0097 3865.5799
42 6827.7963 6276.0097
43 5794.2706 6827.7963
44 7942.4728 5794.2706
45 3552.0242 7942.4728
46 1395.9397 3552.0242
47 -1504.1158 1395.9397
48 475.4229 -1504.1158
49 -2605.4953 475.4229
50 -7625.2404 -2605.4953
51 -7699.0269 -7625.2404
52 -1136.0075 -7699.0269
53 -2625.6339 -1136.0075
54 -3025.3468 -2625.6339
55 280.7475 -3025.3468
56 -7306.2939 280.7475
57 -7054.2695 -7306.2939
58 -8780.5795 -7054.2695
59 -8854.3142 -8780.5795
60 -8290.0955 -8854.3142
61 -9890.6219 -8290.0955
62 -12071.5997 -9890.6219
63 -9393.0571 -12071.5997
64 -3880.4048 -9393.0571
65 1413.1764 -3880.4048
66 1989.8726 1413.1764
67 4318.2094 1989.8726
68 2546.0176 4318.2094
69 2256.1555 2546.0176
70 7546.7131 2256.1555
71 8717.7399 7546.7131
72 11282.0916 8717.7399
73 14168.8779 11282.0916
74 15156.2927 14168.8779
75 NA 15156.2927
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 1107.1534 -1015.7160
[2,] 774.8153 1107.1534
[3,] 9245.3258 774.8153
[4,] 2552.8818 9245.3258
[5,] 4861.9943 2552.8818
[6,] 4567.9963 4861.9943
[7,] 2150.1252 4567.9963
[8,] 4555.7867 2150.1252
[9,] 2993.3220 4555.7867
[10,] 1637.3861 2993.3220
[11,] 1297.0945 1637.3861
[12,] 1806.9008 1297.0945
[13,] 1407.5119 1806.9008
[14,] 5926.8749 1407.5119
[15,] 7443.3713 5926.8749
[16,] 207.2118 7443.3713
[17,] -11896.6661 207.2118
[18,] -12215.5416 -11896.6661
[19,] -10956.1878 -12215.5416
[20,] -6789.0778 -10956.1878
[21,] -1016.6045 -6789.0778
[22,] -3080.4166 -1016.6045
[23,] -1894.3420 -3080.4166
[24,] -4645.6402 -1894.3420
[25,] -4354.0344 -4645.6402
[26,] -1062.8592 -4354.0344
[27,] 1912.6503 -1062.8592
[28,] -1609.2611 1912.6503
[29,] 1971.1196 -1609.2611
[30,] 1855.2232 1971.1196
[31,] -1587.1649 1855.2232
[32,] -948.9054 -1587.1649
[33,] -730.6276 -948.9054
[34,] 1280.9572 -730.6276
[35,] 2237.9376 1280.9572
[36,] 387.0364 2237.9376
[37,] 166.6083 387.0364
[38,] -1098.2838 166.6083
[39,] -1509.2634 -1098.2838
[40,] 3865.5799 -1509.2634
[41,] 6276.0097 3865.5799
[42,] 6827.7963 6276.0097
[43,] 5794.2706 6827.7963
[44,] 7942.4728 5794.2706
[45,] 3552.0242 7942.4728
[46,] 1395.9397 3552.0242
[47,] -1504.1158 1395.9397
[48,] 475.4229 -1504.1158
[49,] -2605.4953 475.4229
[50,] -7625.2404 -2605.4953
[51,] -7699.0269 -7625.2404
[52,] -1136.0075 -7699.0269
[53,] -2625.6339 -1136.0075
[54,] -3025.3468 -2625.6339
[55,] 280.7475 -3025.3468
[56,] -7306.2939 280.7475
[57,] -7054.2695 -7306.2939
[58,] -8780.5795 -7054.2695
[59,] -8854.3142 -8780.5795
[60,] -8290.0955 -8854.3142
[61,] -9890.6219 -8290.0955
[62,] -12071.5997 -9890.6219
[63,] -9393.0571 -12071.5997
[64,] -3880.4048 -9393.0571
[65,] 1413.1764 -3880.4048
[66,] 1989.8726 1413.1764
[67,] 4318.2094 1989.8726
[68,] 2546.0176 4318.2094
[69,] 2256.1555 2546.0176
[70,] 7546.7131 2256.1555
[71,] 8717.7399 7546.7131
[72,] 11282.0916 8717.7399
[73,] 14168.8779 11282.0916
[74,] 15156.2927 14168.8779
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 1107.1534 -1015.7160
2 774.8153 1107.1534
3 9245.3258 774.8153
4 2552.8818 9245.3258
5 4861.9943 2552.8818
6 4567.9963 4861.9943
7 2150.1252 4567.9963
8 4555.7867 2150.1252
9 2993.3220 4555.7867
10 1637.3861 2993.3220
11 1297.0945 1637.3861
12 1806.9008 1297.0945
13 1407.5119 1806.9008
14 5926.8749 1407.5119
15 7443.3713 5926.8749
16 207.2118 7443.3713
17 -11896.6661 207.2118
18 -12215.5416 -11896.6661
19 -10956.1878 -12215.5416
20 -6789.0778 -10956.1878
21 -1016.6045 -6789.0778
22 -3080.4166 -1016.6045
23 -1894.3420 -3080.4166
24 -4645.6402 -1894.3420
25 -4354.0344 -4645.6402
26 -1062.8592 -4354.0344
27 1912.6503 -1062.8592
28 -1609.2611 1912.6503
29 1971.1196 -1609.2611
30 1855.2232 1971.1196
31 -1587.1649 1855.2232
32 -948.9054 -1587.1649
33 -730.6276 -948.9054
34 1280.9572 -730.6276
35 2237.9376 1280.9572
36 387.0364 2237.9376
37 166.6083 387.0364
38 -1098.2838 166.6083
39 -1509.2634 -1098.2838
40 3865.5799 -1509.2634
41 6276.0097 3865.5799
42 6827.7963 6276.0097
43 5794.2706 6827.7963
44 7942.4728 5794.2706
45 3552.0242 7942.4728
46 1395.9397 3552.0242
47 -1504.1158 1395.9397
48 475.4229 -1504.1158
49 -2605.4953 475.4229
50 -7625.2404 -2605.4953
51 -7699.0269 -7625.2404
52 -1136.0075 -7699.0269
53 -2625.6339 -1136.0075
54 -3025.3468 -2625.6339
55 280.7475 -3025.3468
56 -7306.2939 280.7475
57 -7054.2695 -7306.2939
58 -8780.5795 -7054.2695
59 -8854.3142 -8780.5795
60 -8290.0955 -8854.3142
61 -9890.6219 -8290.0955
62 -12071.5997 -9890.6219
63 -9393.0571 -12071.5997
64 -3880.4048 -9393.0571
65 1413.1764 -3880.4048
66 1989.8726 1413.1764
67 4318.2094 1989.8726
68 2546.0176 4318.2094
69 2256.1555 2546.0176
70 7546.7131 2256.1555
71 8717.7399 7546.7131
72 11282.0916 8717.7399
73 14168.8779 11282.0916
74 15156.2927 14168.8779
> 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/70k6n1258904733.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/8edko1258904733.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/9oju81258904733.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/10j73g1258904733.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/11jimi1258904733.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/12wzrt1258904733.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/13ubfn1258904733.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/14zw961258904733.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/15dp701258904733.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/16pmxm1258904733.tab")
+ }
> system("convert tmp/1dson1258904733.ps tmp/1dson1258904733.png")
> system("convert tmp/2d3sm1258904733.ps tmp/2d3sm1258904733.png")
> system("convert tmp/3ymb91258904733.ps tmp/3ymb91258904733.png")
> system("convert tmp/4wc8n1258904733.ps tmp/4wc8n1258904733.png")
> system("convert tmp/59q9s1258904733.ps tmp/59q9s1258904733.png")
> system("convert tmp/68hnx1258904733.ps tmp/68hnx1258904733.png")
> system("convert tmp/70k6n1258904733.ps tmp/70k6n1258904733.png")
> system("convert tmp/8edko1258904733.ps tmp/8edko1258904733.png")
> system("convert tmp/9oju81258904733.ps tmp/9oju81258904733.png")
> system("convert tmp/10j73g1258904733.ps tmp/10j73g1258904733.png")
>
>
> proc.time()
user system elapsed
2.516 1.583 3.551