R version 2.15.2 (2012-10-26) -- "Trick or Treat"
Copyright (C) 2012 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: i686-pc-linux-gnu (32-bit)
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(6.8
+ ,225
+ ,0.442
+ ,0.672
+ ,9.2
+ ,6.3
+ ,180
+ ,0.435
+ ,0.797
+ ,11.7
+ ,6.4
+ ,190
+ ,0.456
+ ,0.761
+ ,15.8
+ ,6.2
+ ,180
+ ,0.416
+ ,0.651
+ ,8.6
+ ,6.9
+ ,205
+ ,0.449
+ ,0.9
+ ,23.2
+ ,6.4
+ ,225
+ ,0.431
+ ,0.78
+ ,27.4
+ ,6.3
+ ,185
+ ,0.487
+ ,0.771
+ ,9.3
+ ,6.8
+ ,235
+ ,0.469
+ ,0.75
+ ,16
+ ,6.9
+ ,235
+ ,0.435
+ ,0.818
+ ,4.7
+ ,6.7
+ ,210
+ ,0.48
+ ,0.825
+ ,12.5
+ ,6.9
+ ,245
+ ,0.516
+ ,0.632
+ ,20.1
+ ,6.9
+ ,245
+ ,0.493
+ ,0.757
+ ,9.1
+ ,6.3
+ ,185
+ ,0.374
+ ,0.709
+ ,8.1
+ ,6.1
+ ,185
+ ,0.424
+ ,0.782
+ ,8.6
+ ,6.2
+ ,180
+ ,0.441
+ ,0.775
+ ,20.3
+ ,6.8
+ ,220
+ ,0.503
+ ,0.88
+ ,25
+ ,6.5
+ ,194
+ ,0.503
+ ,0.833
+ ,19.2
+ ,7.6
+ ,225
+ ,0.425
+ ,0.571
+ ,3.3
+ ,6.3
+ ,210
+ ,0.371
+ ,0.816
+ ,11.2
+ ,7.1
+ ,240
+ ,0.504
+ ,0.714
+ ,10.5
+ ,6.8
+ ,225
+ ,0.4
+ ,0.765
+ ,10.1
+ ,7.3
+ ,263
+ ,0.482
+ ,0.655
+ ,7.2
+ ,6.4
+ ,210
+ ,0.475
+ ,0.244
+ ,13.6
+ ,6.8
+ ,235
+ ,0.428
+ ,0.728
+ ,9
+ ,7.2
+ ,230
+ ,0.559
+ ,0.721
+ ,24.6
+ ,6.4
+ ,190
+ ,0.441
+ ,0.757
+ ,12.6
+ ,6.6
+ ,220
+ ,0.492
+ ,0.747
+ ,5.6
+ ,6.8
+ ,210
+ ,0.402
+ ,0.739
+ ,8.7
+ ,6.1
+ ,180
+ ,0.415
+ ,0.713
+ ,7.7
+ ,6.5
+ ,235
+ ,0.492
+ ,0.742
+ ,24.1
+ ,6.4
+ ,185
+ ,0.484
+ ,0.861
+ ,11.7
+ ,6
+ ,175
+ ,0.387
+ ,0.721
+ ,7.7
+ ,6
+ ,192
+ ,0.436
+ ,0.785
+ ,9.6
+ ,7.3
+ ,263
+ ,0.482
+ ,0.655
+ ,7.2
+ ,6.1
+ ,180
+ ,0.34
+ ,0.821
+ ,12.3
+ ,6.7
+ ,240
+ ,0.516
+ ,0.728
+ ,8.9
+ ,6.4
+ ,210
+ ,0.475
+ ,0.846
+ ,13.6
+ ,5.8
+ ,160
+ ,0.412
+ ,0.813
+ ,11.2
+ ,6.9
+ ,230
+ ,0.411
+ ,0.595
+ ,2.8
+ ,7
+ ,245
+ ,0.407
+ ,0.573
+ ,3.2
+ ,7.3
+ ,228
+ ,0.445
+ ,0.726
+ ,9.4
+ ,5.9
+ ,155
+ ,0.291
+ ,0.707
+ ,11.9
+ ,6.2
+ ,200
+ ,0.449
+ ,0.804
+ ,15.4
+ ,6.8
+ ,235
+ ,0.546
+ ,0.784
+ ,7.4
+ ,7
+ ,235
+ ,0.48
+ ,0.744
+ ,18.9
+ ,5.9
+ ,105
+ ,0.359
+ ,0.839
+ ,7.9
+ ,6.1
+ ,180
+ ,0.528
+ ,0.79
+ ,12.2
+ ,5.7
+ ,185
+ ,0.352
+ ,0.701
+ ,11
+ ,7.1
+ ,245
+ ,0.414
+ ,0.778
+ ,2.8
+ ,5.8
+ ,180
+ ,0.425
+ ,0.872
+ ,11.8
+ ,7.4
+ ,240
+ ,0.599
+ ,0.713
+ ,17.1
+ ,6.8
+ ,225
+ ,0.482
+ ,0.701
+ ,11.6
+ ,6.8
+ ,215
+ ,0.457
+ ,0.734
+ ,5.8
+ ,7
+ ,230
+ ,0.435
+ ,0.764
+ ,8.3)
+ ,dim=c(5
+ ,54)
+ ,dimnames=list(c('X1'
+ ,'X2'
+ ,'X3'
+ ,'X4'
+ ,'X5
')
+ ,1:54))
> y <- array(NA,dim=c(5,54),dimnames=list(c('X1','X2','X3','X4','X5
'),1:54))
> for (i in 1:dim(x)[1])
+ {
+ for (j in 1:dim(x)[2])
+ {
+ y[i,j] <- as.numeric(x[i,j])
+ }
+ }
> par3 = 'No Linear Trend'
> par2 = 'Do not include Seasonal Dummies'
> par1 = '5'
> par3 <- 'No Linear Trend'
> par2 <- 'Do not include Seasonal Dummies'
> par1 <- '5'
> #'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, 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
X5\r X1 X2 X3 X4
1 9.2 6.8 225 0.442 0.672
2 11.7 6.3 180 0.435 0.797
3 15.8 6.4 190 0.456 0.761
4 8.6 6.2 180 0.416 0.651
5 23.2 6.9 205 0.449 0.900
6 27.4 6.4 225 0.431 0.780
7 9.3 6.3 185 0.487 0.771
8 16.0 6.8 235 0.469 0.750
9 4.7 6.9 235 0.435 0.818
10 12.5 6.7 210 0.480 0.825
11 20.1 6.9 245 0.516 0.632
12 9.1 6.9 245 0.493 0.757
13 8.1 6.3 185 0.374 0.709
14 8.6 6.1 185 0.424 0.782
15 20.3 6.2 180 0.441 0.775
16 25.0 6.8 220 0.503 0.880
17 19.2 6.5 194 0.503 0.833
18 3.3 7.6 225 0.425 0.571
19 11.2 6.3 210 0.371 0.816
20 10.5 7.1 240 0.504 0.714
21 10.1 6.8 225 0.400 0.765
22 7.2 7.3 263 0.482 0.655
23 13.6 6.4 210 0.475 0.244
24 9.0 6.8 235 0.428 0.728
25 24.6 7.2 230 0.559 0.721
26 12.6 6.4 190 0.441 0.757
27 5.6 6.6 220 0.492 0.747
28 8.7 6.8 210 0.402 0.739
29 7.7 6.1 180 0.415 0.713
30 24.1 6.5 235 0.492 0.742
31 11.7 6.4 185 0.484 0.861
32 7.7 6.0 175 0.387 0.721
33 9.6 6.0 192 0.436 0.785
34 7.2 7.3 263 0.482 0.655
35 12.3 6.1 180 0.340 0.821
36 8.9 6.7 240 0.516 0.728
37 13.6 6.4 210 0.475 0.846
38 11.2 5.8 160 0.412 0.813
39 2.8 6.9 230 0.411 0.595
40 3.2 7.0 245 0.407 0.573
41 9.4 7.3 228 0.445 0.726
42 11.9 5.9 155 0.291 0.707
43 15.4 6.2 200 0.449 0.804
44 7.4 6.8 235 0.546 0.784
45 18.9 7.0 235 0.480 0.744
46 7.9 5.9 105 0.359 0.839
47 12.2 6.1 180 0.528 0.790
48 11.0 5.7 185 0.352 0.701
49 2.8 7.1 245 0.414 0.778
50 11.8 5.8 180 0.425 0.872
51 17.1 7.4 240 0.599 0.713
52 11.6 6.8 225 0.482 0.701
53 5.8 6.8 215 0.457 0.734
54 8.3 7.0 230 0.435 0.764
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X1 X2 X3 X4
4.148707 -3.690499 0.009458 47.940199 11.371019
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-8.966 -3.545 -1.187 2.613 15.211
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.148707 14.855006 0.279 0.78121
X1 -3.690499 2.970780 -1.242 0.22005
X2 0.009458 0.046297 0.204 0.83897
X3 47.940199 15.709131 3.052 0.00367 **
X4 11.371019 7.868536 1.445 0.15479
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 5.411 on 49 degrees of freedom
Multiple R-squared: 0.2223, Adjusted R-squared: 0.1588
F-statistic: 3.501 on 4 and 49 DF, p-value: 0.01364
> 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.3466745 0.69334893 0.65332554
[2,] 0.9696735 0.06065299 0.03032650
[3,] 0.9423893 0.11522136 0.05761068
[4,] 0.9646046 0.07079089 0.03539545
[5,] 0.9706414 0.05871711 0.02935856
[6,] 0.9485522 0.10289561 0.05144781
[7,] 0.9410528 0.11789434 0.05894717
[8,] 0.9529799 0.09404025 0.04702013
[9,] 0.9734743 0.05305139 0.02652570
[10,] 0.9633986 0.07320281 0.03660140
[11,] 0.9464511 0.10709772 0.05354886
[12,] 0.9225840 0.15483190 0.07741595
[13,] 0.8972987 0.20540255 0.10270128
[14,] 0.8567308 0.28653831 0.14326916
[15,] 0.8186356 0.36272886 0.18136443
[16,] 0.8245053 0.35098937 0.17549469
[17,] 0.7677238 0.46455242 0.23227621
[18,] 0.9007571 0.19848582 0.09924291
[19,] 0.8639124 0.27217518 0.13608759
[20,] 0.9255508 0.14889848 0.07444924
[21,] 0.8900778 0.21984442 0.10992221
[22,] 0.8616793 0.27664134 0.13832067
[23,] 0.9717370 0.05652598 0.02826299
[24,] 0.9613871 0.07722578 0.03861289
[25,] 0.9419608 0.11607834 0.05803917
[26,] 0.9234015 0.15319693 0.07659846
[27,] 0.8943145 0.21137109 0.10568554
[28,] 0.8713855 0.25722900 0.12861450
[29,] 0.8560701 0.28785976 0.14392988
[30,] 0.7990309 0.40193816 0.20096908
[31,] 0.7236130 0.55277395 0.27638698
[32,] 0.7032993 0.59340137 0.29670069
[33,] 0.7462493 0.50750135 0.25375067
[34,] 0.6460079 0.70798429 0.35399214
[35,] 0.5872508 0.82549831 0.41274915
[36,] 0.5520392 0.89592164 0.44796082
[37,] 0.5905766 0.81884688 0.40942344
[38,] 0.8524166 0.29516684 0.14758342
[39,] 0.7989757 0.40204854 0.20102427
> postscript(file="/var/wessaorg/rcomp/tmp/1xzow1353441367.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/wessaorg/rcomp/tmp/2wktn1353441367.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/wessaorg/rcomp/tmp/3kk2q1353441367.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/wessaorg/rcomp/tmp/4ojrf1353441367.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/wessaorg/rcomp/tmp/5ki391353441367.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 = 54
Frequency = 1
1 2 3 4 5 6
-0.81235894 -0.81777389 2.95930395 -1.71579120 10.81768634 15.21071354
7 8 9 10 11 12
-5.46231003 3.71173160 -6.36248103 -1.30102549 7.17478785 -4.14396498
13 14 15 16 17 18
-0.54006433 -4.00525851 7.37569743 9.74540920 3.61861728 -1.79650334
19 20 21 22 23 24
1.25059576 -1.99696124 1.04362464 -3.05083143 5.53808797 -1.07255781
25 26 27 28 29 30
9.85036516 0.52389102 -8.35300287 -0.01473239 -3.64190410 9.69292545
31 32 33 34 35 36
-3.57283126 -2.71230430 -4.04991307 -3.05083143 3.32554075 -5.80763753
37 38 39 40 41 42
-1.30726563 -2.05316600 -4.52888666 -3.45979040 0.44665959 6.06926834
43 44 45 46 47 48
1.57324712 -8.96627839 6.89071535 -2.21871685 -5.43471509 1.29128875
49 50 51 52 53 54
-6.15738083 -2.93644788 1.16724059 -0.65972646 -5.54188054 -1.73210379
> postscript(file="/var/wessaorg/rcomp/tmp/6bznr1353441367.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 = 54
Frequency = 1
lag(myerror, k = 1) myerror
0 -0.81235894 NA
1 -0.81777389 -0.81235894
2 2.95930395 -0.81777389
3 -1.71579120 2.95930395
4 10.81768634 -1.71579120
5 15.21071354 10.81768634
6 -5.46231003 15.21071354
7 3.71173160 -5.46231003
8 -6.36248103 3.71173160
9 -1.30102549 -6.36248103
10 7.17478785 -1.30102549
11 -4.14396498 7.17478785
12 -0.54006433 -4.14396498
13 -4.00525851 -0.54006433
14 7.37569743 -4.00525851
15 9.74540920 7.37569743
16 3.61861728 9.74540920
17 -1.79650334 3.61861728
18 1.25059576 -1.79650334
19 -1.99696124 1.25059576
20 1.04362464 -1.99696124
21 -3.05083143 1.04362464
22 5.53808797 -3.05083143
23 -1.07255781 5.53808797
24 9.85036516 -1.07255781
25 0.52389102 9.85036516
26 -8.35300287 0.52389102
27 -0.01473239 -8.35300287
28 -3.64190410 -0.01473239
29 9.69292545 -3.64190410
30 -3.57283126 9.69292545
31 -2.71230430 -3.57283126
32 -4.04991307 -2.71230430
33 -3.05083143 -4.04991307
34 3.32554075 -3.05083143
35 -5.80763753 3.32554075
36 -1.30726563 -5.80763753
37 -2.05316600 -1.30726563
38 -4.52888666 -2.05316600
39 -3.45979040 -4.52888666
40 0.44665959 -3.45979040
41 6.06926834 0.44665959
42 1.57324712 6.06926834
43 -8.96627839 1.57324712
44 6.89071535 -8.96627839
45 -2.21871685 6.89071535
46 -5.43471509 -2.21871685
47 1.29128875 -5.43471509
48 -6.15738083 1.29128875
49 -2.93644788 -6.15738083
50 1.16724059 -2.93644788
51 -0.65972646 1.16724059
52 -5.54188054 -0.65972646
53 -1.73210379 -5.54188054
54 NA -1.73210379
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -0.81777389 -0.81235894
[2,] 2.95930395 -0.81777389
[3,] -1.71579120 2.95930395
[4,] 10.81768634 -1.71579120
[5,] 15.21071354 10.81768634
[6,] -5.46231003 15.21071354
[7,] 3.71173160 -5.46231003
[8,] -6.36248103 3.71173160
[9,] -1.30102549 -6.36248103
[10,] 7.17478785 -1.30102549
[11,] -4.14396498 7.17478785
[12,] -0.54006433 -4.14396498
[13,] -4.00525851 -0.54006433
[14,] 7.37569743 -4.00525851
[15,] 9.74540920 7.37569743
[16,] 3.61861728 9.74540920
[17,] -1.79650334 3.61861728
[18,] 1.25059576 -1.79650334
[19,] -1.99696124 1.25059576
[20,] 1.04362464 -1.99696124
[21,] -3.05083143 1.04362464
[22,] 5.53808797 -3.05083143
[23,] -1.07255781 5.53808797
[24,] 9.85036516 -1.07255781
[25,] 0.52389102 9.85036516
[26,] -8.35300287 0.52389102
[27,] -0.01473239 -8.35300287
[28,] -3.64190410 -0.01473239
[29,] 9.69292545 -3.64190410
[30,] -3.57283126 9.69292545
[31,] -2.71230430 -3.57283126
[32,] -4.04991307 -2.71230430
[33,] -3.05083143 -4.04991307
[34,] 3.32554075 -3.05083143
[35,] -5.80763753 3.32554075
[36,] -1.30726563 -5.80763753
[37,] -2.05316600 -1.30726563
[38,] -4.52888666 -2.05316600
[39,] -3.45979040 -4.52888666
[40,] 0.44665959 -3.45979040
[41,] 6.06926834 0.44665959
[42,] 1.57324712 6.06926834
[43,] -8.96627839 1.57324712
[44,] 6.89071535 -8.96627839
[45,] -2.21871685 6.89071535
[46,] -5.43471509 -2.21871685
[47,] 1.29128875 -5.43471509
[48,] -6.15738083 1.29128875
[49,] -2.93644788 -6.15738083
[50,] 1.16724059 -2.93644788
[51,] -0.65972646 1.16724059
[52,] -5.54188054 -0.65972646
[53,] -1.73210379 -5.54188054
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -0.81777389 -0.81235894
2 2.95930395 -0.81777389
3 -1.71579120 2.95930395
4 10.81768634 -1.71579120
5 15.21071354 10.81768634
6 -5.46231003 15.21071354
7 3.71173160 -5.46231003
8 -6.36248103 3.71173160
9 -1.30102549 -6.36248103
10 7.17478785 -1.30102549
11 -4.14396498 7.17478785
12 -0.54006433 -4.14396498
13 -4.00525851 -0.54006433
14 7.37569743 -4.00525851
15 9.74540920 7.37569743
16 3.61861728 9.74540920
17 -1.79650334 3.61861728
18 1.25059576 -1.79650334
19 -1.99696124 1.25059576
20 1.04362464 -1.99696124
21 -3.05083143 1.04362464
22 5.53808797 -3.05083143
23 -1.07255781 5.53808797
24 9.85036516 -1.07255781
25 0.52389102 9.85036516
26 -8.35300287 0.52389102
27 -0.01473239 -8.35300287
28 -3.64190410 -0.01473239
29 9.69292545 -3.64190410
30 -3.57283126 9.69292545
31 -2.71230430 -3.57283126
32 -4.04991307 -2.71230430
33 -3.05083143 -4.04991307
34 3.32554075 -3.05083143
35 -5.80763753 3.32554075
36 -1.30726563 -5.80763753
37 -2.05316600 -1.30726563
38 -4.52888666 -2.05316600
39 -3.45979040 -4.52888666
40 0.44665959 -3.45979040
41 6.06926834 0.44665959
42 1.57324712 6.06926834
43 -8.96627839 1.57324712
44 6.89071535 -8.96627839
45 -2.21871685 6.89071535
46 -5.43471509 -2.21871685
47 1.29128875 -5.43471509
48 -6.15738083 1.29128875
49 -2.93644788 -6.15738083
50 1.16724059 -2.93644788
51 -0.65972646 1.16724059
52 -5.54188054 -0.65972646
53 -1.73210379 -5.54188054
> 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/wessaorg/rcomp/tmp/7bkz81353441367.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/wessaorg/rcomp/tmp/8tnqs1353441367.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/wessaorg/rcomp/tmp/91pas1353441367.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/wessaorg/rcomp/tmp/10ygaz1353441367.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/wessaorg/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/wessaorg/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/wessaorg/rcomp/tmp/11ia5t1353441367.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/wessaorg/rcomp/tmp/12gczs1353441367.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/wessaorg/rcomp/tmp/13wlt01353441368.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/wessaorg/rcomp/tmp/14w0pl1353441368.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/wessaorg/rcomp/tmp/159xkw1353441368.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/wessaorg/rcomp/tmp/16uky91353441368.tab")
+ }
>
> try(system("convert tmp/1xzow1353441367.ps tmp/1xzow1353441367.png",intern=TRUE))
character(0)
> try(system("convert tmp/2wktn1353441367.ps tmp/2wktn1353441367.png",intern=TRUE))
character(0)
> try(system("convert tmp/3kk2q1353441367.ps tmp/3kk2q1353441367.png",intern=TRUE))
character(0)
> try(system("convert tmp/4ojrf1353441367.ps tmp/4ojrf1353441367.png",intern=TRUE))
character(0)
> try(system("convert tmp/5ki391353441367.ps tmp/5ki391353441367.png",intern=TRUE))
character(0)
> try(system("convert tmp/6bznr1353441367.ps tmp/6bznr1353441367.png",intern=TRUE))
character(0)
> try(system("convert tmp/7bkz81353441367.ps tmp/7bkz81353441367.png",intern=TRUE))
character(0)
> try(system("convert tmp/8tnqs1353441367.ps tmp/8tnqs1353441367.png",intern=TRUE))
character(0)
> try(system("convert tmp/91pas1353441367.ps tmp/91pas1353441367.png",intern=TRUE))
character(0)
> try(system("convert tmp/10ygaz1353441367.ps tmp/10ygaz1353441367.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
6.326 1.181 7.574