R version 2.7.0 (2008-04-22)
Copyright (C) 2008 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(101.3,11554.5,102,13182.1,109.2,14800.1,88.6,12150.7,94.3,14478.2,98.3,13253.9,86.4,12036.8,80.6,12653.2,104.1,14035.4,108.2,14571.4,93.4,15400.9,71.9,14283.2,94.1,14485.3,94.9,14196.3,96.4,15559.1,91.1,13767.4,84.4,14634,86.4,14381.1,88,12509.9,75.1,12122.3,109.7,13122.3,103,13908.7,82.1,13456.5,68,12441.6,96.4,12953,94.3,13057.2,90,14350.1,88,13830.2,76.1,13755.5,82.5,13574.4,81.4,12802.6,66.5,11737.3,97.2,13850.2,94.1,15081.8,80.7,13653.3,70.5,14019.1,87.8,13962,89.5,13768.7,99.6,14747.1,84.2,13858.1,75.1,13188,92,13693.1,80.8,12970,73.1,11392.8,99.8,13985.2,90,14994.7,83.1,13584.7,72.4,14257.8,78.8,13553.4,87.3,14007.3,91,16535.8,80.1,14721.4,73.6,13664.6,86.4,16405.9,74.5,13829.4,71.2,13735.6,92.4,15870.5,81.5,15962.4,85.3,15744.1,69.9,16083.7,84.2,14863.9,90.7,15533.1,100.3,17473.1,79.4,15925.5,84.8,15573.7,92.9,17495,81.6,14155.8,76,14913.9,98.7,17250.4,89.1,15879.8,88.7,17647.8,67.1,17749.9),dim=c(2,72),dimnames=list(c('textiel','Invoer'),1:72))
> y <- array(NA,dim=c(2,72),dimnames=list(c('textiel','Invoer'),1:72))
> 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 = '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
textiel Invoer M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 101.3 11554.5 1 0 0 0 0 0 0 0 0 0 0
2 102.0 13182.1 0 1 0 0 0 0 0 0 0 0 0
3 109.2 14800.1 0 0 1 0 0 0 0 0 0 0 0
4 88.6 12150.7 0 0 0 1 0 0 0 0 0 0 0
5 94.3 14478.2 0 0 0 0 1 0 0 0 0 0 0
6 98.3 13253.9 0 0 0 0 0 1 0 0 0 0 0
7 86.4 12036.8 0 0 0 0 0 0 1 0 0 0 0
8 80.6 12653.2 0 0 0 0 0 0 0 1 0 0 0
9 104.1 14035.4 0 0 0 0 0 0 0 0 1 0 0
10 108.2 14571.4 0 0 0 0 0 0 0 0 0 1 0
11 93.4 15400.9 0 0 0 0 0 0 0 0 0 0 1
12 71.9 14283.2 0 0 0 0 0 0 0 0 0 0 0
13 94.1 14485.3 1 0 0 0 0 0 0 0 0 0 0
14 94.9 14196.3 0 1 0 0 0 0 0 0 0 0 0
15 96.4 15559.1 0 0 1 0 0 0 0 0 0 0 0
16 91.1 13767.4 0 0 0 1 0 0 0 0 0 0 0
17 84.4 14634.0 0 0 0 0 1 0 0 0 0 0 0
18 86.4 14381.1 0 0 0 0 0 1 0 0 0 0 0
19 88.0 12509.9 0 0 0 0 0 0 1 0 0 0 0
20 75.1 12122.3 0 0 0 0 0 0 0 1 0 0 0
21 109.7 13122.3 0 0 0 0 0 0 0 0 1 0 0
22 103.0 13908.7 0 0 0 0 0 0 0 0 0 1 0
23 82.1 13456.5 0 0 0 0 0 0 0 0 0 0 1
24 68.0 12441.6 0 0 0 0 0 0 0 0 0 0 0
25 96.4 12953.0 1 0 0 0 0 0 0 0 0 0 0
26 94.3 13057.2 0 1 0 0 0 0 0 0 0 0 0
27 90.0 14350.1 0 0 1 0 0 0 0 0 0 0 0
28 88.0 13830.2 0 0 0 1 0 0 0 0 0 0 0
29 76.1 13755.5 0 0 0 0 1 0 0 0 0 0 0
30 82.5 13574.4 0 0 0 0 0 1 0 0 0 0 0
31 81.4 12802.6 0 0 0 0 0 0 1 0 0 0 0
32 66.5 11737.3 0 0 0 0 0 0 0 1 0 0 0
33 97.2 13850.2 0 0 0 0 0 0 0 0 1 0 0
34 94.1 15081.8 0 0 0 0 0 0 0 0 0 1 0
35 80.7 13653.3 0 0 0 0 0 0 0 0 0 0 1
36 70.5 14019.1 0 0 0 0 0 0 0 0 0 0 0
37 87.8 13962.0 1 0 0 0 0 0 0 0 0 0 0
38 89.5 13768.7 0 1 0 0 0 0 0 0 0 0 0
39 99.6 14747.1 0 0 1 0 0 0 0 0 0 0 0
40 84.2 13858.1 0 0 0 1 0 0 0 0 0 0 0
41 75.1 13188.0 0 0 0 0 1 0 0 0 0 0 0
42 92.0 13693.1 0 0 0 0 0 1 0 0 0 0 0
43 80.8 12970.0 0 0 0 0 0 0 1 0 0 0 0
44 73.1 11392.8 0 0 0 0 0 0 0 1 0 0 0
45 99.8 13985.2 0 0 0 0 0 0 0 0 1 0 0
46 90.0 14994.7 0 0 0 0 0 0 0 0 0 1 0
47 83.1 13584.7 0 0 0 0 0 0 0 0 0 0 1
48 72.4 14257.8 0 0 0 0 0 0 0 0 0 0 0
49 78.8 13553.4 1 0 0 0 0 0 0 0 0 0 0
50 87.3 14007.3 0 1 0 0 0 0 0 0 0 0 0
51 91.0 16535.8 0 0 1 0 0 0 0 0 0 0 0
52 80.1 14721.4 0 0 0 1 0 0 0 0 0 0 0
53 73.6 13664.6 0 0 0 0 1 0 0 0 0 0 0
54 86.4 16405.9 0 0 0 0 0 1 0 0 0 0 0
55 74.5 13829.4 0 0 0 0 0 0 1 0 0 0 0
56 71.2 13735.6 0 0 0 0 0 0 0 1 0 0 0
57 92.4 15870.5 0 0 0 0 0 0 0 0 1 0 0
58 81.5 15962.4 0 0 0 0 0 0 0 0 0 1 0
59 85.3 15744.1 0 0 0 0 0 0 0 0 0 0 1
60 69.9 16083.7 0 0 0 0 0 0 0 0 0 0 0
61 84.2 14863.9 1 0 0 0 0 0 0 0 0 0 0
62 90.7 15533.1 0 1 0 0 0 0 0 0 0 0 0
63 100.3 17473.1 0 0 1 0 0 0 0 0 0 0 0
64 79.4 15925.5 0 0 0 1 0 0 0 0 0 0 0
65 84.8 15573.7 0 0 0 0 1 0 0 0 0 0 0
66 92.9 17495.0 0 0 0 0 0 1 0 0 0 0 0
67 81.6 14155.8 0 0 0 0 0 0 1 0 0 0 0
68 76.0 14913.9 0 0 0 0 0 0 0 1 0 0 0
69 98.7 17250.4 0 0 0 0 0 0 0 0 1 0 0
70 89.1 15879.8 0 0 0 0 0 0 0 0 0 1 0
71 88.7 17647.8 0 0 0 0 0 0 0 0 0 0 1
72 67.1 17749.9 0 0 0 0 0 0 0 0 0 0 0
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Invoer M1 M2 M3 M4
83.1762254 -0.0008922 19.3569099 22.3930421 28.4718013 14.5853362
M5 M6 M7 M8 M9 M10
10.8900853 19.7785899 10.5840997 1.9573022 30.2427447 24.5824880
M11
15.6802839
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-12.0173 -4.2100 -0.1858 4.1896 13.4416
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 83.1762254 9.3745667 8.873 1.86e-12 ***
Invoer -0.0008922 0.0006097 -1.463 0.148695
M1 19.3569099 3.6545822 5.297 1.83e-06 ***
M2 22.3930421 3.6122570 6.199 6.03e-08 ***
M3 28.4718013 3.6058518 7.896 8.20e-11 ***
M4 14.5853362 3.6052186 4.046 0.000154 ***
M5 10.8900853 3.5930910 3.031 0.003618 **
M6 19.7785899 3.5750264 5.532 7.60e-07 ***
M7 10.5840997 3.7317510 2.836 0.006245 **
M8 1.9573022 3.7865607 0.517 0.607153
M9 30.2427447 3.5757762 8.458 9.25e-12 ***
M10 24.5824880 3.5785536 6.869 4.50e-09 ***
M11 15.6802839 3.5756388 4.385 4.84e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 6.192 on 59 degrees of freedom
Multiple R-squared: 0.7144, Adjusted R-squared: 0.6563
F-statistic: 12.3 on 12 and 59 DF, p-value: 5.261e-12
> 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.6626664 0.67466713 0.337333564
[2,] 0.7187410 0.56251807 0.281259036
[3,] 0.7144418 0.57111634 0.285558169
[4,] 0.6430822 0.71383566 0.356917831
[5,] 0.6197528 0.76049445 0.380247223
[6,] 0.6134568 0.77308636 0.386543182
[7,] 0.7453746 0.50925070 0.254625350
[8,] 0.8954563 0.20908747 0.104543735
[9,] 0.8672985 0.26540297 0.132701483
[10,] 0.9095126 0.18097472 0.090487359
[11,] 0.9023113 0.19537743 0.097688716
[12,] 0.9627136 0.07457276 0.037286381
[13,] 0.9603716 0.07925684 0.039628420
[14,] 0.9789614 0.04207711 0.021038553
[15,] 0.9870421 0.02591574 0.012957868
[16,] 0.9832648 0.03347038 0.016735189
[17,] 0.9902395 0.01952090 0.009760452
[18,] 0.9895693 0.02086146 0.010430729
[19,] 0.9949887 0.01002261 0.005011305
[20,] 0.9942434 0.01151321 0.005756604
[21,] 0.9897342 0.02053157 0.010265785
[22,] 0.9923725 0.01525501 0.007627504
[23,] 0.9891852 0.02162961 0.010814804
[24,] 0.9851948 0.02961047 0.014805237
[25,] 0.9825058 0.03498847 0.017494235
[26,] 0.9807988 0.03840237 0.019201187
[27,] 0.9736066 0.05278685 0.026393424
[28,] 0.9632079 0.07358427 0.036792133
[29,] 0.9419500 0.11609993 0.058049966
[30,] 0.9470288 0.10594243 0.052971214
[31,] 0.9580306 0.08393888 0.041969439
[32,] 0.9311729 0.13765421 0.068827107
[33,] 0.9725096 0.05498076 0.027490380
[34,] 0.9694828 0.06103437 0.030517185
[35,] 0.9488015 0.10239694 0.051198468
[36,] 0.9448300 0.11033991 0.055169956
[37,] 0.9205175 0.15896493 0.079482465
[38,] 0.9098510 0.18029791 0.090148956
[39,] 0.8608417 0.27831653 0.139158264
[40,] 0.8477661 0.30446789 0.152233945
[41,] 0.7345442 0.53091166 0.265455828
> postscript(file="/var/www/html/rcomp/tmp/19yis1229763789.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/2qlt31229763789.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/3aks11229763789.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/4oqz71229763789.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/55dc11229763789.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 = 72
Frequency = 1
1 2 3 4 5 6
9.0755942 8.1915792 10.7563723 1.6790875 13.1508945 7.1700902
7 8 9 10 11 12
3.3787043 6.7554435 3.2031764 13.4416432 8.2839132 1.4670041
13 14 15 16 17 18
4.4904044 1.9964313 -1.3664608 5.6214799 3.3898966 -3.7242411
19 20 21 22 23 24
5.4007961 0.7817835 7.9885240 7.6503935 -4.7508476 -4.0760403
25 26 27 28 29 30
5.4233123 0.3801455 -8.8451101 2.5775090 -5.6938863 -8.3439652
31 32 33 34 35 36
-0.9380619 -8.1617070 -3.8620559 -0.2029865 -5.9752659 -0.1686215
37 38 39 40 41 42
-2.2764750 -3.7850662 1.1090866 -1.1975991 -7.2002002 1.2619370
43 44 45 46 47 48
-1.3887105 -1.8690641 -1.1416112 -4.3806957 -3.6364697 1.9443427
49 50 51 52 53 54
-11.6410210 -5.7721913 -5.8950656 -4.5273774 -8.2749857 -1.9177487
55 56 57 58 59 60
-6.9219683 -1.6788576 -6.8595784 -12.0173301 0.4901105 1.0733798
61 62 63 64 65 66
-5.0718150 -1.0108984 4.2411777 -4.1530998 4.6282810 5.5539279
67 68 69 70 71 72
0.4692403 4.1724018 0.6715450 -4.4910244 5.5885594 -0.2400648
> postscript(file="/var/www/html/rcomp/tmp/60c0b1229763789.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 = 72
Frequency = 1
lag(myerror, k = 1) myerror
0 9.0755942 NA
1 8.1915792 9.0755942
2 10.7563723 8.1915792
3 1.6790875 10.7563723
4 13.1508945 1.6790875
5 7.1700902 13.1508945
6 3.3787043 7.1700902
7 6.7554435 3.3787043
8 3.2031764 6.7554435
9 13.4416432 3.2031764
10 8.2839132 13.4416432
11 1.4670041 8.2839132
12 4.4904044 1.4670041
13 1.9964313 4.4904044
14 -1.3664608 1.9964313
15 5.6214799 -1.3664608
16 3.3898966 5.6214799
17 -3.7242411 3.3898966
18 5.4007961 -3.7242411
19 0.7817835 5.4007961
20 7.9885240 0.7817835
21 7.6503935 7.9885240
22 -4.7508476 7.6503935
23 -4.0760403 -4.7508476
24 5.4233123 -4.0760403
25 0.3801455 5.4233123
26 -8.8451101 0.3801455
27 2.5775090 -8.8451101
28 -5.6938863 2.5775090
29 -8.3439652 -5.6938863
30 -0.9380619 -8.3439652
31 -8.1617070 -0.9380619
32 -3.8620559 -8.1617070
33 -0.2029865 -3.8620559
34 -5.9752659 -0.2029865
35 -0.1686215 -5.9752659
36 -2.2764750 -0.1686215
37 -3.7850662 -2.2764750
38 1.1090866 -3.7850662
39 -1.1975991 1.1090866
40 -7.2002002 -1.1975991
41 1.2619370 -7.2002002
42 -1.3887105 1.2619370
43 -1.8690641 -1.3887105
44 -1.1416112 -1.8690641
45 -4.3806957 -1.1416112
46 -3.6364697 -4.3806957
47 1.9443427 -3.6364697
48 -11.6410210 1.9443427
49 -5.7721913 -11.6410210
50 -5.8950656 -5.7721913
51 -4.5273774 -5.8950656
52 -8.2749857 -4.5273774
53 -1.9177487 -8.2749857
54 -6.9219683 -1.9177487
55 -1.6788576 -6.9219683
56 -6.8595784 -1.6788576
57 -12.0173301 -6.8595784
58 0.4901105 -12.0173301
59 1.0733798 0.4901105
60 -5.0718150 1.0733798
61 -1.0108984 -5.0718150
62 4.2411777 -1.0108984
63 -4.1530998 4.2411777
64 4.6282810 -4.1530998
65 5.5539279 4.6282810
66 0.4692403 5.5539279
67 4.1724018 0.4692403
68 0.6715450 4.1724018
69 -4.4910244 0.6715450
70 5.5885594 -4.4910244
71 -0.2400648 5.5885594
72 NA -0.2400648
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 8.1915792 9.0755942
[2,] 10.7563723 8.1915792
[3,] 1.6790875 10.7563723
[4,] 13.1508945 1.6790875
[5,] 7.1700902 13.1508945
[6,] 3.3787043 7.1700902
[7,] 6.7554435 3.3787043
[8,] 3.2031764 6.7554435
[9,] 13.4416432 3.2031764
[10,] 8.2839132 13.4416432
[11,] 1.4670041 8.2839132
[12,] 4.4904044 1.4670041
[13,] 1.9964313 4.4904044
[14,] -1.3664608 1.9964313
[15,] 5.6214799 -1.3664608
[16,] 3.3898966 5.6214799
[17,] -3.7242411 3.3898966
[18,] 5.4007961 -3.7242411
[19,] 0.7817835 5.4007961
[20,] 7.9885240 0.7817835
[21,] 7.6503935 7.9885240
[22,] -4.7508476 7.6503935
[23,] -4.0760403 -4.7508476
[24,] 5.4233123 -4.0760403
[25,] 0.3801455 5.4233123
[26,] -8.8451101 0.3801455
[27,] 2.5775090 -8.8451101
[28,] -5.6938863 2.5775090
[29,] -8.3439652 -5.6938863
[30,] -0.9380619 -8.3439652
[31,] -8.1617070 -0.9380619
[32,] -3.8620559 -8.1617070
[33,] -0.2029865 -3.8620559
[34,] -5.9752659 -0.2029865
[35,] -0.1686215 -5.9752659
[36,] -2.2764750 -0.1686215
[37,] -3.7850662 -2.2764750
[38,] 1.1090866 -3.7850662
[39,] -1.1975991 1.1090866
[40,] -7.2002002 -1.1975991
[41,] 1.2619370 -7.2002002
[42,] -1.3887105 1.2619370
[43,] -1.8690641 -1.3887105
[44,] -1.1416112 -1.8690641
[45,] -4.3806957 -1.1416112
[46,] -3.6364697 -4.3806957
[47,] 1.9443427 -3.6364697
[48,] -11.6410210 1.9443427
[49,] -5.7721913 -11.6410210
[50,] -5.8950656 -5.7721913
[51,] -4.5273774 -5.8950656
[52,] -8.2749857 -4.5273774
[53,] -1.9177487 -8.2749857
[54,] -6.9219683 -1.9177487
[55,] -1.6788576 -6.9219683
[56,] -6.8595784 -1.6788576
[57,] -12.0173301 -6.8595784
[58,] 0.4901105 -12.0173301
[59,] 1.0733798 0.4901105
[60,] -5.0718150 1.0733798
[61,] -1.0108984 -5.0718150
[62,] 4.2411777 -1.0108984
[63,] -4.1530998 4.2411777
[64,] 4.6282810 -4.1530998
[65,] 5.5539279 4.6282810
[66,] 0.4692403 5.5539279
[67,] 4.1724018 0.4692403
[68,] 0.6715450 4.1724018
[69,] -4.4910244 0.6715450
[70,] 5.5885594 -4.4910244
[71,] -0.2400648 5.5885594
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 8.1915792 9.0755942
2 10.7563723 8.1915792
3 1.6790875 10.7563723
4 13.1508945 1.6790875
5 7.1700902 13.1508945
6 3.3787043 7.1700902
7 6.7554435 3.3787043
8 3.2031764 6.7554435
9 13.4416432 3.2031764
10 8.2839132 13.4416432
11 1.4670041 8.2839132
12 4.4904044 1.4670041
13 1.9964313 4.4904044
14 -1.3664608 1.9964313
15 5.6214799 -1.3664608
16 3.3898966 5.6214799
17 -3.7242411 3.3898966
18 5.4007961 -3.7242411
19 0.7817835 5.4007961
20 7.9885240 0.7817835
21 7.6503935 7.9885240
22 -4.7508476 7.6503935
23 -4.0760403 -4.7508476
24 5.4233123 -4.0760403
25 0.3801455 5.4233123
26 -8.8451101 0.3801455
27 2.5775090 -8.8451101
28 -5.6938863 2.5775090
29 -8.3439652 -5.6938863
30 -0.9380619 -8.3439652
31 -8.1617070 -0.9380619
32 -3.8620559 -8.1617070
33 -0.2029865 -3.8620559
34 -5.9752659 -0.2029865
35 -0.1686215 -5.9752659
36 -2.2764750 -0.1686215
37 -3.7850662 -2.2764750
38 1.1090866 -3.7850662
39 -1.1975991 1.1090866
40 -7.2002002 -1.1975991
41 1.2619370 -7.2002002
42 -1.3887105 1.2619370
43 -1.8690641 -1.3887105
44 -1.1416112 -1.8690641
45 -4.3806957 -1.1416112
46 -3.6364697 -4.3806957
47 1.9443427 -3.6364697
48 -11.6410210 1.9443427
49 -5.7721913 -11.6410210
50 -5.8950656 -5.7721913
51 -4.5273774 -5.8950656
52 -8.2749857 -4.5273774
53 -1.9177487 -8.2749857
54 -6.9219683 -1.9177487
55 -1.6788576 -6.9219683
56 -6.8595784 -1.6788576
57 -12.0173301 -6.8595784
58 0.4901105 -12.0173301
59 1.0733798 0.4901105
60 -5.0718150 1.0733798
61 -1.0108984 -5.0718150
62 4.2411777 -1.0108984
63 -4.1530998 4.2411777
64 4.6282810 -4.1530998
65 5.5539279 4.6282810
66 0.4692403 5.5539279
67 4.1724018 0.4692403
68 0.6715450 4.1724018
69 -4.4910244 0.6715450
70 5.5885594 -4.4910244
71 -0.2400648 5.5885594
> 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/7w9wo1229763789.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/8fzoz1229763789.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/9613y1229763789.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/10jtwy1229763789.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/11qodk1229763789.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/12u8he1229763789.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/13jluv1229763789.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/149rbx1229763789.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/15xcpf1229763789.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/16xf051229763789.tab")
+ }
>
> system("convert tmp/19yis1229763789.ps tmp/19yis1229763789.png")
> system("convert tmp/2qlt31229763789.ps tmp/2qlt31229763789.png")
> system("convert tmp/3aks11229763789.ps tmp/3aks11229763789.png")
> system("convert tmp/4oqz71229763789.ps tmp/4oqz71229763789.png")
> system("convert tmp/55dc11229763789.ps tmp/55dc11229763789.png")
> system("convert tmp/60c0b1229763789.ps tmp/60c0b1229763789.png")
> system("convert tmp/7w9wo1229763789.ps tmp/7w9wo1229763789.png")
> system("convert tmp/8fzoz1229763789.ps tmp/8fzoz1229763789.png")
> system("convert tmp/9613y1229763789.ps tmp/9613y1229763789.png")
> system("convert tmp/10jtwy1229763789.ps tmp/10jtwy1229763789.png")
>
>
> proc.time()
user system elapsed
5.252 2.761 5.616