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(11554.5
+ ,7.5
+ ,13182.1
+ ,7.2
+ ,14800.1
+ ,6.9
+ ,12150.7
+ ,6.7
+ ,14478.2
+ ,6.4
+ ,13253.9
+ ,6.3
+ ,12036.8
+ ,6.8
+ ,12653.2
+ ,7.3
+ ,14035.4
+ ,7.1
+ ,14571.4
+ ,7.1
+ ,15400.9
+ ,6.8
+ ,14283.2
+ ,6.5
+ ,14485.3
+ ,6.3
+ ,14196.3
+ ,6.1
+ ,15559.1
+ ,6.1
+ ,13767.4
+ ,6.3
+ ,14634
+ ,6.3
+ ,14381.1
+ ,6
+ ,12509.9
+ ,6.2
+ ,12122.3
+ ,6.4
+ ,13122.3
+ ,6.8
+ ,13908.7
+ ,7.5
+ ,13456.5
+ ,7.5
+ ,12441.6
+ ,7.6
+ ,12953
+ ,7.6
+ ,13057.2
+ ,7.4
+ ,14350.1
+ ,7.3
+ ,13830.2
+ ,7.1
+ ,13755.5
+ ,6.9
+ ,13574.4
+ ,6.8
+ ,12802.6
+ ,7.5
+ ,11737.3
+ ,7.6
+ ,13850.2
+ ,7.8
+ ,15081.8
+ ,8
+ ,13653.3
+ ,8.1
+ ,14019.1
+ ,8.2
+ ,13962
+ ,8.3
+ ,13768.7
+ ,8.2
+ ,14747.1
+ ,8
+ ,13858.1
+ ,7.9
+ ,13188
+ ,7.6
+ ,13693.1
+ ,7.6
+ ,12970
+ ,8.2
+ ,11392.8
+ ,8.3
+ ,13985.2
+ ,8.4
+ ,14994.7
+ ,8.4
+ ,13584.7
+ ,8.4
+ ,14257.8
+ ,8.6
+ ,13553.4
+ ,8.9
+ ,14007.3
+ ,8.8
+ ,16535.8
+ ,8.3
+ ,14721.4
+ ,7.5
+ ,13664.6
+ ,7.2
+ ,16805.9
+ ,7.5
+ ,13829.4
+ ,8.8
+ ,13735.6
+ ,9.3
+ ,15870.5
+ ,9.3
+ ,15962.4
+ ,8.7
+ ,15744.1
+ ,8.2
+ ,16083.7
+ ,8.3
+ ,14863.9
+ ,8.5
+ ,15533.1
+ ,8.6
+ ,17473.1
+ ,8.6
+ ,15925.5
+ ,8.2
+ ,15573.7
+ ,8.1
+ ,17495
+ ,8
+ ,14155.8
+ ,8.6
+ ,14913.9
+ ,8.7
+ ,17250.4
+ ,8.8
+ ,15879.8
+ ,8.5
+ ,17647.8
+ ,8.4
+ ,17749.9
+ ,8.5
+ ,17111.8
+ ,8.7
+ ,16934.8
+ ,8.7
+ ,20280
+ ,8.6
+ ,16238.2
+ ,8.5
+ ,17896.1
+ ,8.3
+ ,18089.3
+ ,8.1
+ ,15660
+ ,8.2
+ ,16162.4
+ ,8.1
+ ,17850.1
+ ,8.1
+ ,18520.4
+ ,7.9
+ ,18524.7
+ ,7.9
+ ,16843.7
+ ,7.9)
+ ,dim=c(2
+ ,84)
+ ,dimnames=list(c('Invoer'
+ ,'Werkloosheid')
+ ,1:84))
> y <- array(NA,dim=c(2,84),dimnames=list(c('Invoer','Werkloosheid'),1:84))
> 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 = '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
Invoer Werkloosheid
1 11554.5 7.5
2 13182.1 7.2
3 14800.1 6.9
4 12150.7 6.7
5 14478.2 6.4
6 13253.9 6.3
7 12036.8 6.8
8 12653.2 7.3
9 14035.4 7.1
10 14571.4 7.1
11 15400.9 6.8
12 14283.2 6.5
13 14485.3 6.3
14 14196.3 6.1
15 15559.1 6.1
16 13767.4 6.3
17 14634.0 6.3
18 14381.1 6.0
19 12509.9 6.2
20 12122.3 6.4
21 13122.3 6.8
22 13908.7 7.5
23 13456.5 7.5
24 12441.6 7.6
25 12953.0 7.6
26 13057.2 7.4
27 14350.1 7.3
28 13830.2 7.1
29 13755.5 6.9
30 13574.4 6.8
31 12802.6 7.5
32 11737.3 7.6
33 13850.2 7.8
34 15081.8 8.0
35 13653.3 8.1
36 14019.1 8.2
37 13962.0 8.3
38 13768.7 8.2
39 14747.1 8.0
40 13858.1 7.9
41 13188.0 7.6
42 13693.1 7.6
43 12970.0 8.2
44 11392.8 8.3
45 13985.2 8.4
46 14994.7 8.4
47 13584.7 8.4
48 14257.8 8.6
49 13553.4 8.9
50 14007.3 8.8
51 16535.8 8.3
52 14721.4 7.5
53 13664.6 7.2
54 16805.9 7.5
55 13829.4 8.8
56 13735.6 9.3
57 15870.5 9.3
58 15962.4 8.7
59 15744.1 8.2
60 16083.7 8.3
61 14863.9 8.5
62 15533.1 8.6
63 17473.1 8.6
64 15925.5 8.2
65 15573.7 8.1
66 17495.0 8.0
67 14155.8 8.6
68 14913.9 8.7
69 17250.4 8.8
70 15879.8 8.5
71 17647.8 8.4
72 17749.9 8.5
73 17111.8 8.7
74 16934.8 8.7
75 20280.0 8.6
76 16238.2 8.5
77 17896.1 8.3
78 18089.3 8.1
79 15660.0 8.2
80 16162.4 8.1
81 17850.1 8.1
82 18520.4 7.9
83 18524.7 7.9
84 16843.7 7.9
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Werkloosheid
7986.8 872.6
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-3836.7 -1332.6 -225.5 1027.6 4788.7
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 7986.8 1729.2 4.619 1.41e-05 ***
Werkloosheid 872.6 221.8 3.935 0.000174 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1686 on 82 degrees of freedom
Multiple R-squared: 0.1588, Adjusted R-squared: 0.1485
F-statistic: 15.48 on 1 and 82 DF, p-value: 0.0001737
> 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,] 3.906936e-01 0.7813871704 0.6093064
[2,] 2.843810e-01 0.5687619248 0.7156190
[3,] 2.376312e-01 0.4752623369 0.7623688
[4,] 1.460917e-01 0.2921834316 0.8539083
[5,] 1.245970e-01 0.2491940317 0.8754030
[6,] 1.328293e-01 0.2656585956 0.8671707
[7,] 1.760896e-01 0.3521792551 0.8239104
[8,] 1.175782e-01 0.2351564333 0.8824218
[9,] 7.543444e-02 0.1508688876 0.9245656
[10,] 4.666058e-02 0.0933211660 0.9533394
[11,] 3.962707e-02 0.0792541305 0.9603729
[12,] 2.540823e-02 0.0508164600 0.9745918
[13,] 1.558424e-02 0.0311684704 0.9844158
[14,] 9.660278e-03 0.0193205552 0.9903397
[15,] 1.419757e-02 0.0283951442 0.9858024
[16,] 1.966704e-02 0.0393340745 0.9803330
[17,] 1.217467e-02 0.0243493469 0.9878253
[18,] 8.978865e-03 0.0179577302 0.9910211
[19,] 5.495241e-03 0.0109904819 0.9945048
[20,] 3.947032e-03 0.0078940648 0.9960530
[21,] 2.442098e-03 0.0048841963 0.9975579
[22,] 1.450956e-03 0.0029019117 0.9985490
[23,] 1.149896e-03 0.0022997914 0.9988501
[24,] 6.502229e-04 0.0013004459 0.9993498
[25,] 3.440698e-04 0.0006881396 0.9996559
[26,] 1.786420e-04 0.0003572840 0.9998214
[27,] 1.178330e-04 0.0002356661 0.9998822
[28,] 2.032202e-04 0.0004064403 0.9997968
[29,] 1.724908e-04 0.0003449817 0.9998275
[30,] 3.965433e-04 0.0007930866 0.9996035
[31,] 2.900157e-04 0.0005800313 0.9997100
[32,] 2.261062e-04 0.0004522124 0.9997739
[33,] 1.684363e-04 0.0003368727 0.9998316
[34,] 1.196898e-04 0.0002393796 0.9998803
[35,] 1.099627e-04 0.0002199254 0.9998900
[36,] 7.693861e-05 0.0001538772 0.9999231
[37,] 7.605421e-05 0.0001521084 0.9999239
[38,] 6.717179e-05 0.0001343436 0.9999328
[39,] 8.646434e-05 0.0001729287 0.9999135
[40,] 1.076829e-03 0.0021536588 0.9989232
[41,] 1.107270e-03 0.0022145393 0.9988927
[42,] 1.377327e-03 0.0027546540 0.9986227
[43,] 1.663962e-03 0.0033279237 0.9983360
[44,] 1.668853e-03 0.0033377052 0.9983311
[45,] 2.031050e-03 0.0040621001 0.9979689
[46,] 2.357081e-03 0.0047141629 0.9976429
[47,] 7.000478e-03 0.0140009568 0.9929995
[48,] 8.831770e-03 0.0176635402 0.9911682
[49,] 4.257061e-02 0.0851412198 0.9574294
[50,] 9.432608e-02 0.1886521520 0.9056739
[51,] 1.142773e-01 0.2285546893 0.8857227
[52,] 1.210262e-01 0.2420523348 0.8789738
[53,] 1.229583e-01 0.2459165729 0.8770417
[54,] 1.246430e-01 0.2492859485 0.8753570
[55,] 1.359035e-01 0.2718070876 0.8640965
[56,] 1.418619e-01 0.2837237821 0.8581381
[57,] 1.628270e-01 0.3256539216 0.8371730
[58,] 1.600097e-01 0.3200194537 0.8399903
[59,] 2.081541e-01 0.4163081672 0.7918459
[60,] 2.102795e-01 0.4205589958 0.7897205
[61,] 2.426813e-01 0.4853625844 0.7573187
[62,] 2.679281e-01 0.5358562988 0.7320719
[63,] 4.383574e-01 0.8767148920 0.5616426
[64,] 5.580467e-01 0.8839066517 0.4419533
[65,] 5.265842e-01 0.9468315542 0.4734158
[66,] 5.614658e-01 0.8770684250 0.4385342
[67,] 5.346921e-01 0.9306157905 0.4653079
[68,] 4.985239e-01 0.9970478634 0.5014761
[69,] 4.340469e-01 0.8680938222 0.5659531
[70,] 3.841925e-01 0.7683850587 0.6158075
[71,] 8.688278e-01 0.2623444641 0.1311722
[72,] 7.880147e-01 0.4239705236 0.2119853
[73,] 8.262108e-01 0.3475783740 0.1737892
[74,] 8.343148e-01 0.3313704484 0.1656852
[75,] 7.228488e-01 0.5543024559 0.2771512
> postscript(file="/var/www/html/rcomp/tmp/1mbmy1228659414.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/2r2io1228659414.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/3yrqb1228659414.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/4a6th1228659414.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/514kv1228659414.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 = 84
Frequency = 1
1 2 3 4 5 6
-2976.926326 -1087.541109 792.244107 -1682.632415 906.652801 -230.385460
7 8 9 10 11 12
-1883.794154 -1703.702848 -146.979371 389.020629 1480.305846 624.391062
13 14 15 16 17 18
1001.014540 886.538018 2249.338018 283.114540 1149.714540 1158.599757
19 20 21 22 23 24
-887.123721 -1449.247199 -798.294154 -622.726326 -1074.926326 -2177.088065
25 26 27 28 29 30
-1665.688065 -1386.964587 -6.802848 -352.179371 -252.355893 -346.194154
31 32 33 34 35 36
-1728.826326 -2881.388065 -943.011542 114.064980 -1401.696759 -1123.158498
37 38 39 40 41 42
-1267.520236 -1373.558498 -220.635020 -1022.373281 -1430.688065 -925.588065
43 44 45 46 47 48
-2172.258498 -3836.720236 -1331.581975 -322.081975 -1732.081975 -1233.505453
49 50 51 52 53 54
-2199.690669 -1658.528931 1306.279764 189.973674 -605.041109 2274.473674
55 56 57 58 59 60
-1836.428931 -2366.537625 -231.637625 383.832808 601.841502 854.179764
61 62 63 64 65 66
-540.143714 41.794547 1981.794547 783.241502 518.703241 2527.264980
67 68 69 70 71 72
-1335.505453 -664.667192 1584.571069 475.756286 2331.018025 2345.856286
73 74 75 76 77 78
1533.232808 1356.232808 4788.694547 834.156286 2666.579764 3034.303241
79 80 81 82 83 84
517.741502 1107.403241 2795.103241 3639.926719 3644.226719 1963.226719
> postscript(file="/var/www/html/rcomp/tmp/6rnak1228659414.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 = 84
Frequency = 1
lag(myerror, k = 1) myerror
0 -2976.926326 NA
1 -1087.541109 -2976.926326
2 792.244107 -1087.541109
3 -1682.632415 792.244107
4 906.652801 -1682.632415
5 -230.385460 906.652801
6 -1883.794154 -230.385460
7 -1703.702848 -1883.794154
8 -146.979371 -1703.702848
9 389.020629 -146.979371
10 1480.305846 389.020629
11 624.391062 1480.305846
12 1001.014540 624.391062
13 886.538018 1001.014540
14 2249.338018 886.538018
15 283.114540 2249.338018
16 1149.714540 283.114540
17 1158.599757 1149.714540
18 -887.123721 1158.599757
19 -1449.247199 -887.123721
20 -798.294154 -1449.247199
21 -622.726326 -798.294154
22 -1074.926326 -622.726326
23 -2177.088065 -1074.926326
24 -1665.688065 -2177.088065
25 -1386.964587 -1665.688065
26 -6.802848 -1386.964587
27 -352.179371 -6.802848
28 -252.355893 -352.179371
29 -346.194154 -252.355893
30 -1728.826326 -346.194154
31 -2881.388065 -1728.826326
32 -943.011542 -2881.388065
33 114.064980 -943.011542
34 -1401.696759 114.064980
35 -1123.158498 -1401.696759
36 -1267.520236 -1123.158498
37 -1373.558498 -1267.520236
38 -220.635020 -1373.558498
39 -1022.373281 -220.635020
40 -1430.688065 -1022.373281
41 -925.588065 -1430.688065
42 -2172.258498 -925.588065
43 -3836.720236 -2172.258498
44 -1331.581975 -3836.720236
45 -322.081975 -1331.581975
46 -1732.081975 -322.081975
47 -1233.505453 -1732.081975
48 -2199.690669 -1233.505453
49 -1658.528931 -2199.690669
50 1306.279764 -1658.528931
51 189.973674 1306.279764
52 -605.041109 189.973674
53 2274.473674 -605.041109
54 -1836.428931 2274.473674
55 -2366.537625 -1836.428931
56 -231.637625 -2366.537625
57 383.832808 -231.637625
58 601.841502 383.832808
59 854.179764 601.841502
60 -540.143714 854.179764
61 41.794547 -540.143714
62 1981.794547 41.794547
63 783.241502 1981.794547
64 518.703241 783.241502
65 2527.264980 518.703241
66 -1335.505453 2527.264980
67 -664.667192 -1335.505453
68 1584.571069 -664.667192
69 475.756286 1584.571069
70 2331.018025 475.756286
71 2345.856286 2331.018025
72 1533.232808 2345.856286
73 1356.232808 1533.232808
74 4788.694547 1356.232808
75 834.156286 4788.694547
76 2666.579764 834.156286
77 3034.303241 2666.579764
78 517.741502 3034.303241
79 1107.403241 517.741502
80 2795.103241 1107.403241
81 3639.926719 2795.103241
82 3644.226719 3639.926719
83 1963.226719 3644.226719
84 NA 1963.226719
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -1087.541109 -2976.926326
[2,] 792.244107 -1087.541109
[3,] -1682.632415 792.244107
[4,] 906.652801 -1682.632415
[5,] -230.385460 906.652801
[6,] -1883.794154 -230.385460
[7,] -1703.702848 -1883.794154
[8,] -146.979371 -1703.702848
[9,] 389.020629 -146.979371
[10,] 1480.305846 389.020629
[11,] 624.391062 1480.305846
[12,] 1001.014540 624.391062
[13,] 886.538018 1001.014540
[14,] 2249.338018 886.538018
[15,] 283.114540 2249.338018
[16,] 1149.714540 283.114540
[17,] 1158.599757 1149.714540
[18,] -887.123721 1158.599757
[19,] -1449.247199 -887.123721
[20,] -798.294154 -1449.247199
[21,] -622.726326 -798.294154
[22,] -1074.926326 -622.726326
[23,] -2177.088065 -1074.926326
[24,] -1665.688065 -2177.088065
[25,] -1386.964587 -1665.688065
[26,] -6.802848 -1386.964587
[27,] -352.179371 -6.802848
[28,] -252.355893 -352.179371
[29,] -346.194154 -252.355893
[30,] -1728.826326 -346.194154
[31,] -2881.388065 -1728.826326
[32,] -943.011542 -2881.388065
[33,] 114.064980 -943.011542
[34,] -1401.696759 114.064980
[35,] -1123.158498 -1401.696759
[36,] -1267.520236 -1123.158498
[37,] -1373.558498 -1267.520236
[38,] -220.635020 -1373.558498
[39,] -1022.373281 -220.635020
[40,] -1430.688065 -1022.373281
[41,] -925.588065 -1430.688065
[42,] -2172.258498 -925.588065
[43,] -3836.720236 -2172.258498
[44,] -1331.581975 -3836.720236
[45,] -322.081975 -1331.581975
[46,] -1732.081975 -322.081975
[47,] -1233.505453 -1732.081975
[48,] -2199.690669 -1233.505453
[49,] -1658.528931 -2199.690669
[50,] 1306.279764 -1658.528931
[51,] 189.973674 1306.279764
[52,] -605.041109 189.973674
[53,] 2274.473674 -605.041109
[54,] -1836.428931 2274.473674
[55,] -2366.537625 -1836.428931
[56,] -231.637625 -2366.537625
[57,] 383.832808 -231.637625
[58,] 601.841502 383.832808
[59,] 854.179764 601.841502
[60,] -540.143714 854.179764
[61,] 41.794547 -540.143714
[62,] 1981.794547 41.794547
[63,] 783.241502 1981.794547
[64,] 518.703241 783.241502
[65,] 2527.264980 518.703241
[66,] -1335.505453 2527.264980
[67,] -664.667192 -1335.505453
[68,] 1584.571069 -664.667192
[69,] 475.756286 1584.571069
[70,] 2331.018025 475.756286
[71,] 2345.856286 2331.018025
[72,] 1533.232808 2345.856286
[73,] 1356.232808 1533.232808
[74,] 4788.694547 1356.232808
[75,] 834.156286 4788.694547
[76,] 2666.579764 834.156286
[77,] 3034.303241 2666.579764
[78,] 517.741502 3034.303241
[79,] 1107.403241 517.741502
[80,] 2795.103241 1107.403241
[81,] 3639.926719 2795.103241
[82,] 3644.226719 3639.926719
[83,] 1963.226719 3644.226719
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -1087.541109 -2976.926326
2 792.244107 -1087.541109
3 -1682.632415 792.244107
4 906.652801 -1682.632415
5 -230.385460 906.652801
6 -1883.794154 -230.385460
7 -1703.702848 -1883.794154
8 -146.979371 -1703.702848
9 389.020629 -146.979371
10 1480.305846 389.020629
11 624.391062 1480.305846
12 1001.014540 624.391062
13 886.538018 1001.014540
14 2249.338018 886.538018
15 283.114540 2249.338018
16 1149.714540 283.114540
17 1158.599757 1149.714540
18 -887.123721 1158.599757
19 -1449.247199 -887.123721
20 -798.294154 -1449.247199
21 -622.726326 -798.294154
22 -1074.926326 -622.726326
23 -2177.088065 -1074.926326
24 -1665.688065 -2177.088065
25 -1386.964587 -1665.688065
26 -6.802848 -1386.964587
27 -352.179371 -6.802848
28 -252.355893 -352.179371
29 -346.194154 -252.355893
30 -1728.826326 -346.194154
31 -2881.388065 -1728.826326
32 -943.011542 -2881.388065
33 114.064980 -943.011542
34 -1401.696759 114.064980
35 -1123.158498 -1401.696759
36 -1267.520236 -1123.158498
37 -1373.558498 -1267.520236
38 -220.635020 -1373.558498
39 -1022.373281 -220.635020
40 -1430.688065 -1022.373281
41 -925.588065 -1430.688065
42 -2172.258498 -925.588065
43 -3836.720236 -2172.258498
44 -1331.581975 -3836.720236
45 -322.081975 -1331.581975
46 -1732.081975 -322.081975
47 -1233.505453 -1732.081975
48 -2199.690669 -1233.505453
49 -1658.528931 -2199.690669
50 1306.279764 -1658.528931
51 189.973674 1306.279764
52 -605.041109 189.973674
53 2274.473674 -605.041109
54 -1836.428931 2274.473674
55 -2366.537625 -1836.428931
56 -231.637625 -2366.537625
57 383.832808 -231.637625
58 601.841502 383.832808
59 854.179764 601.841502
60 -540.143714 854.179764
61 41.794547 -540.143714
62 1981.794547 41.794547
63 783.241502 1981.794547
64 518.703241 783.241502
65 2527.264980 518.703241
66 -1335.505453 2527.264980
67 -664.667192 -1335.505453
68 1584.571069 -664.667192
69 475.756286 1584.571069
70 2331.018025 475.756286
71 2345.856286 2331.018025
72 1533.232808 2345.856286
73 1356.232808 1533.232808
74 4788.694547 1356.232808
75 834.156286 4788.694547
76 2666.579764 834.156286
77 3034.303241 2666.579764
78 517.741502 3034.303241
79 1107.403241 517.741502
80 2795.103241 1107.403241
81 3639.926719 2795.103241
82 3644.226719 3639.926719
83 1963.226719 3644.226719
> 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/7sdyb1228659414.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/8r3b61228659414.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/9v1ne1228659414.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/104zwq1228659414.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/11n53d1228659414.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/122tl21228659414.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/13ozwu1228659414.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/14x6hx1228659414.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/15ese41228659414.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/16wvxe1228659414.tab")
+ }
>
> system("convert tmp/1mbmy1228659414.ps tmp/1mbmy1228659414.png")
> system("convert tmp/2r2io1228659414.ps tmp/2r2io1228659414.png")
> system("convert tmp/3yrqb1228659414.ps tmp/3yrqb1228659414.png")
> system("convert tmp/4a6th1228659414.ps tmp/4a6th1228659414.png")
> system("convert tmp/514kv1228659414.ps tmp/514kv1228659414.png")
> system("convert tmp/6rnak1228659414.ps tmp/6rnak1228659414.png")
> system("convert tmp/7sdyb1228659414.ps tmp/7sdyb1228659414.png")
> system("convert tmp/8r3b61228659414.ps tmp/8r3b61228659414.png")
> system("convert tmp/9v1ne1228659414.ps tmp/9v1ne1228659414.png")
> system("convert tmp/104zwq1228659414.ps tmp/104zwq1228659414.png")
>
>
> proc.time()
user system elapsed
5.445 2.742 5.853