R version 2.6.2 (2008-02-08)
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(56421,53152,53536,52408,41454,38271,35306,26414,31917,38030,27534,18387,50556,43901,48572,43899,37532,40357,35489,29027,34485,42598,30306,26451,47460,50104,61465,53726,39477,43895,31481,29896,33842,39120,33702,25094,51442,45594,52518,48564,41745,49585,32747,33379,35645,37034,35681,20972,58552,54955,65540,51570,51145,46641,35704,33253,35193,41668,34865,21210,56126,49231,59723,48103,47472,50497,40059,34149,36860,46356,36577),dim=c(1,71),dimnames=list(c('Registraties'),1:71))
> y <- array(NA,dim=c(1,71),dimnames=list(c('Registraties'),1:71))
> 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
> 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
Registraties M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 56421 1 0 0 0 0 0 0 0 0 0 0
2 53152 0 1 0 0 0 0 0 0 0 0 0
3 53536 0 0 1 0 0 0 0 0 0 0 0
4 52408 0 0 0 1 0 0 0 0 0 0 0
5 41454 0 0 0 0 1 0 0 0 0 0 0
6 38271 0 0 0 0 0 1 0 0 0 0 0
7 35306 0 0 0 0 0 0 1 0 0 0 0
8 26414 0 0 0 0 0 0 0 1 0 0 0
9 31917 0 0 0 0 0 0 0 0 1 0 0
10 38030 0 0 0 0 0 0 0 0 0 1 0
11 27534 0 0 0 0 0 0 0 0 0 0 1
12 18387 0 0 0 0 0 0 0 0 0 0 0
13 50556 1 0 0 0 0 0 0 0 0 0 0
14 43901 0 1 0 0 0 0 0 0 0 0 0
15 48572 0 0 1 0 0 0 0 0 0 0 0
16 43899 0 0 0 1 0 0 0 0 0 0 0
17 37532 0 0 0 0 1 0 0 0 0 0 0
18 40357 0 0 0 0 0 1 0 0 0 0 0
19 35489 0 0 0 0 0 0 1 0 0 0 0
20 29027 0 0 0 0 0 0 0 1 0 0 0
21 34485 0 0 0 0 0 0 0 0 1 0 0
22 42598 0 0 0 0 0 0 0 0 0 1 0
23 30306 0 0 0 0 0 0 0 0 0 0 1
24 26451 0 0 0 0 0 0 0 0 0 0 0
25 47460 1 0 0 0 0 0 0 0 0 0 0
26 50104 0 1 0 0 0 0 0 0 0 0 0
27 61465 0 0 1 0 0 0 0 0 0 0 0
28 53726 0 0 0 1 0 0 0 0 0 0 0
29 39477 0 0 0 0 1 0 0 0 0 0 0
30 43895 0 0 0 0 0 1 0 0 0 0 0
31 31481 0 0 0 0 0 0 1 0 0 0 0
32 29896 0 0 0 0 0 0 0 1 0 0 0
33 33842 0 0 0 0 0 0 0 0 1 0 0
34 39120 0 0 0 0 0 0 0 0 0 1 0
35 33702 0 0 0 0 0 0 0 0 0 0 1
36 25094 0 0 0 0 0 0 0 0 0 0 0
37 51442 1 0 0 0 0 0 0 0 0 0 0
38 45594 0 1 0 0 0 0 0 0 0 0 0
39 52518 0 0 1 0 0 0 0 0 0 0 0
40 48564 0 0 0 1 0 0 0 0 0 0 0
41 41745 0 0 0 0 1 0 0 0 0 0 0
42 49585 0 0 0 0 0 1 0 0 0 0 0
43 32747 0 0 0 0 0 0 1 0 0 0 0
44 33379 0 0 0 0 0 0 0 1 0 0 0
45 35645 0 0 0 0 0 0 0 0 1 0 0
46 37034 0 0 0 0 0 0 0 0 0 1 0
47 35681 0 0 0 0 0 0 0 0 0 0 1
48 20972 0 0 0 0 0 0 0 0 0 0 0
49 58552 1 0 0 0 0 0 0 0 0 0 0
50 54955 0 1 0 0 0 0 0 0 0 0 0
51 65540 0 0 1 0 0 0 0 0 0 0 0
52 51570 0 0 0 1 0 0 0 0 0 0 0
53 51145 0 0 0 0 1 0 0 0 0 0 0
54 46641 0 0 0 0 0 1 0 0 0 0 0
55 35704 0 0 0 0 0 0 1 0 0 0 0
56 33253 0 0 0 0 0 0 0 1 0 0 0
57 35193 0 0 0 0 0 0 0 0 1 0 0
58 41668 0 0 0 0 0 0 0 0 0 1 0
59 34865 0 0 0 0 0 0 0 0 0 0 1
60 21210 0 0 0 0 0 0 0 0 0 0 0
61 56126 1 0 0 0 0 0 0 0 0 0 0
62 49231 0 1 0 0 0 0 0 0 0 0 0
63 59723 0 0 1 0 0 0 0 0 0 0 0
64 48103 0 0 0 1 0 0 0 0 0 0 0
65 47472 0 0 0 0 1 0 0 0 0 0 0
66 50497 0 0 0 0 0 1 0 0 0 0 0
67 40059 0 0 0 0 0 0 1 0 0 0 0
68 34149 0 0 0 0 0 0 0 1 0 0 0
69 36860 0 0 0 0 0 0 0 0 1 0 0
70 46356 0 0 0 0 0 0 0 0 0 1 0
71 36577 0 0 0 0 0 0 0 0 0 0 1
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) M1 M2 M3 M4 M5
22423 31003 27067 34470 27289 20715
M6 M7 M8 M9 M10 M11
22452 12708 8597 12234 18378 10688
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-8320 -2788 175 2698 8648
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 22423 1815 12.357 < 2e-16 ***
M1 31003 2457 12.619 < 2e-16 ***
M2 27067 2457 11.016 6.05e-16 ***
M3 34470 2457 14.029 < 2e-16 ***
M4 27289 2457 11.107 4.35e-16 ***
M5 20715 2457 8.431 1.03e-11 ***
M6 22452 2457 9.138 6.73e-13 ***
M7 12708 2457 5.172 2.90e-06 ***
M8 8597 2457 3.499 0.000895 ***
M9 12234 2457 4.979 5.87e-06 ***
M10 18378 2457 7.480 4.16e-10 ***
M11 10688 2457 4.350 5.47e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 4058 on 59 degrees of freedom
Multiple R-squared: 0.8713, Adjusted R-squared: 0.8473
F-statistic: 36.31 on 11 and 59 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.8167001 0.3665997 0.18329987
[2,] 0.8831188 0.2337624 0.11688118
[3,] 0.8615209 0.2769582 0.13847909
[4,] 0.8272708 0.3454583 0.17272915
[5,] 0.7407984 0.5184033 0.25920164
[6,] 0.6743577 0.6512845 0.32564225
[7,] 0.5913478 0.8173044 0.40865221
[8,] 0.5516047 0.8967906 0.44839529
[9,] 0.5036964 0.9926072 0.49630358
[10,] 0.5987242 0.8025517 0.40127584
[11,] 0.7021015 0.5957970 0.29789849
[12,] 0.6249743 0.7500513 0.37502566
[13,] 0.7940635 0.4118730 0.20593649
[14,] 0.8030881 0.3938239 0.19691195
[15,] 0.8121347 0.3757306 0.18786531
[16,] 0.8188832 0.3622335 0.18111676
[17,] 0.8148912 0.3702177 0.18510883
[18,] 0.7839688 0.4320625 0.21603124
[19,] 0.7262055 0.5475889 0.27379447
[20,] 0.6688903 0.6622194 0.33110972
[21,] 0.6415231 0.7169538 0.35847688
[22,] 0.6076020 0.7847960 0.39239799
[23,] 0.6033028 0.7933944 0.39669719
[24,] 0.6472344 0.7055311 0.35276556
[25,] 0.8342044 0.3315912 0.16579560
[26,] 0.7821528 0.4356943 0.21784716
[27,] 0.8766430 0.2467141 0.12335704
[28,] 0.8980044 0.2039911 0.10199556
[29,] 0.9122240 0.1755521 0.08777604
[30,] 0.8871183 0.2257633 0.11288165
[31,] 0.8396952 0.3206095 0.16030476
[32,] 0.9251122 0.1497757 0.07488784
[33,] 0.8976066 0.2047867 0.10239335
[34,] 0.8470310 0.3059381 0.15296904
[35,] 0.8332198 0.3335604 0.16678021
[36,] 0.8843900 0.2312200 0.11561000
[37,] 0.9496889 0.1006221 0.05031106
[38,] 0.9349422 0.1301156 0.06505780
[39,] 0.9416498 0.1167003 0.05835017
[40,] 0.9315741 0.1368518 0.06842588
[41,] 0.9369698 0.1260604 0.06303022
[42,] 0.8525701 0.2948599 0.14742993
> postscript(file="/var/www/html/rcomp/tmp/1hurs1209843625.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/2xyeb1209843625.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/3zvk61209843625.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/42ylf1209843625.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/5la6f1209843625.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 = 71
Frequency = 1
1 2 3 4 5 6 7
2994.8333 3662.5000 -3356.3333 2696.3333 -1683.5000 -6603.3333 175.0000
8 9 10 11 12 13 14
-4605.6667 -2740.0000 -2771.0000 -5576.8333 -4035.8000 -2870.1667 -5588.5000
15 16 17 18 19 20 21
-8320.3333 -5812.6667 -5605.5000 -4517.3333 358.0000 -1992.6667 -172.0000
22 23 24 25 26 27 28
1797.0000 -2804.8333 4028.2000 -5966.1667 614.5000 4572.6667 4014.3333
29 30 31 32 33 34 35
-3660.5000 -979.3333 -3650.0000 -1123.6667 -815.0000 -1681.0000 591.1667
36 37 38 39 40 41 42
2671.2000 -1984.1667 -3895.5000 -4374.3333 -1147.6667 -1392.5000 4710.6667
43 44 45 46 47 48 49
-2384.0000 2359.3333 988.0000 -3767.0000 2570.1667 -1450.8000 5125.8333
50 51 52 53 54 55 56
5465.5000 8647.6667 1858.3333 8007.5000 1766.6667 573.0000 2233.3333
57 58 59 60 61 62 63
536.0000 867.0000 1754.1667 -1212.8000 2699.8333 -258.5000 2830.6667
64 65 66 67 68 69 70
-1608.6667 4334.5000 5622.6667 4928.0000 3129.3333 2203.0000 5555.0000
71
3466.1667
> postscript(file="/var/www/html/rcomp/tmp/6y7yg1209843625.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 = 71
Frequency = 1
lag(myerror, k = 1) myerror
0 2994.8333 NA
1 3662.5000 2994.8333
2 -3356.3333 3662.5000
3 2696.3333 -3356.3333
4 -1683.5000 2696.3333
5 -6603.3333 -1683.5000
6 175.0000 -6603.3333
7 -4605.6667 175.0000
8 -2740.0000 -4605.6667
9 -2771.0000 -2740.0000
10 -5576.8333 -2771.0000
11 -4035.8000 -5576.8333
12 -2870.1667 -4035.8000
13 -5588.5000 -2870.1667
14 -8320.3333 -5588.5000
15 -5812.6667 -8320.3333
16 -5605.5000 -5812.6667
17 -4517.3333 -5605.5000
18 358.0000 -4517.3333
19 -1992.6667 358.0000
20 -172.0000 -1992.6667
21 1797.0000 -172.0000
22 -2804.8333 1797.0000
23 4028.2000 -2804.8333
24 -5966.1667 4028.2000
25 614.5000 -5966.1667
26 4572.6667 614.5000
27 4014.3333 4572.6667
28 -3660.5000 4014.3333
29 -979.3333 -3660.5000
30 -3650.0000 -979.3333
31 -1123.6667 -3650.0000
32 -815.0000 -1123.6667
33 -1681.0000 -815.0000
34 591.1667 -1681.0000
35 2671.2000 591.1667
36 -1984.1667 2671.2000
37 -3895.5000 -1984.1667
38 -4374.3333 -3895.5000
39 -1147.6667 -4374.3333
40 -1392.5000 -1147.6667
41 4710.6667 -1392.5000
42 -2384.0000 4710.6667
43 2359.3333 -2384.0000
44 988.0000 2359.3333
45 -3767.0000 988.0000
46 2570.1667 -3767.0000
47 -1450.8000 2570.1667
48 5125.8333 -1450.8000
49 5465.5000 5125.8333
50 8647.6667 5465.5000
51 1858.3333 8647.6667
52 8007.5000 1858.3333
53 1766.6667 8007.5000
54 573.0000 1766.6667
55 2233.3333 573.0000
56 536.0000 2233.3333
57 867.0000 536.0000
58 1754.1667 867.0000
59 -1212.8000 1754.1667
60 2699.8333 -1212.8000
61 -258.5000 2699.8333
62 2830.6667 -258.5000
63 -1608.6667 2830.6667
64 4334.5000 -1608.6667
65 5622.6667 4334.5000
66 4928.0000 5622.6667
67 3129.3333 4928.0000
68 2203.0000 3129.3333
69 5555.0000 2203.0000
70 3466.1667 5555.0000
71 NA 3466.1667
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 3662.5000 2994.8333
[2,] -3356.3333 3662.5000
[3,] 2696.3333 -3356.3333
[4,] -1683.5000 2696.3333
[5,] -6603.3333 -1683.5000
[6,] 175.0000 -6603.3333
[7,] -4605.6667 175.0000
[8,] -2740.0000 -4605.6667
[9,] -2771.0000 -2740.0000
[10,] -5576.8333 -2771.0000
[11,] -4035.8000 -5576.8333
[12,] -2870.1667 -4035.8000
[13,] -5588.5000 -2870.1667
[14,] -8320.3333 -5588.5000
[15,] -5812.6667 -8320.3333
[16,] -5605.5000 -5812.6667
[17,] -4517.3333 -5605.5000
[18,] 358.0000 -4517.3333
[19,] -1992.6667 358.0000
[20,] -172.0000 -1992.6667
[21,] 1797.0000 -172.0000
[22,] -2804.8333 1797.0000
[23,] 4028.2000 -2804.8333
[24,] -5966.1667 4028.2000
[25,] 614.5000 -5966.1667
[26,] 4572.6667 614.5000
[27,] 4014.3333 4572.6667
[28,] -3660.5000 4014.3333
[29,] -979.3333 -3660.5000
[30,] -3650.0000 -979.3333
[31,] -1123.6667 -3650.0000
[32,] -815.0000 -1123.6667
[33,] -1681.0000 -815.0000
[34,] 591.1667 -1681.0000
[35,] 2671.2000 591.1667
[36,] -1984.1667 2671.2000
[37,] -3895.5000 -1984.1667
[38,] -4374.3333 -3895.5000
[39,] -1147.6667 -4374.3333
[40,] -1392.5000 -1147.6667
[41,] 4710.6667 -1392.5000
[42,] -2384.0000 4710.6667
[43,] 2359.3333 -2384.0000
[44,] 988.0000 2359.3333
[45,] -3767.0000 988.0000
[46,] 2570.1667 -3767.0000
[47,] -1450.8000 2570.1667
[48,] 5125.8333 -1450.8000
[49,] 5465.5000 5125.8333
[50,] 8647.6667 5465.5000
[51,] 1858.3333 8647.6667
[52,] 8007.5000 1858.3333
[53,] 1766.6667 8007.5000
[54,] 573.0000 1766.6667
[55,] 2233.3333 573.0000
[56,] 536.0000 2233.3333
[57,] 867.0000 536.0000
[58,] 1754.1667 867.0000
[59,] -1212.8000 1754.1667
[60,] 2699.8333 -1212.8000
[61,] -258.5000 2699.8333
[62,] 2830.6667 -258.5000
[63,] -1608.6667 2830.6667
[64,] 4334.5000 -1608.6667
[65,] 5622.6667 4334.5000
[66,] 4928.0000 5622.6667
[67,] 3129.3333 4928.0000
[68,] 2203.0000 3129.3333
[69,] 5555.0000 2203.0000
[70,] 3466.1667 5555.0000
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 3662.5000 2994.8333
2 -3356.3333 3662.5000
3 2696.3333 -3356.3333
4 -1683.5000 2696.3333
5 -6603.3333 -1683.5000
6 175.0000 -6603.3333
7 -4605.6667 175.0000
8 -2740.0000 -4605.6667
9 -2771.0000 -2740.0000
10 -5576.8333 -2771.0000
11 -4035.8000 -5576.8333
12 -2870.1667 -4035.8000
13 -5588.5000 -2870.1667
14 -8320.3333 -5588.5000
15 -5812.6667 -8320.3333
16 -5605.5000 -5812.6667
17 -4517.3333 -5605.5000
18 358.0000 -4517.3333
19 -1992.6667 358.0000
20 -172.0000 -1992.6667
21 1797.0000 -172.0000
22 -2804.8333 1797.0000
23 4028.2000 -2804.8333
24 -5966.1667 4028.2000
25 614.5000 -5966.1667
26 4572.6667 614.5000
27 4014.3333 4572.6667
28 -3660.5000 4014.3333
29 -979.3333 -3660.5000
30 -3650.0000 -979.3333
31 -1123.6667 -3650.0000
32 -815.0000 -1123.6667
33 -1681.0000 -815.0000
34 591.1667 -1681.0000
35 2671.2000 591.1667
36 -1984.1667 2671.2000
37 -3895.5000 -1984.1667
38 -4374.3333 -3895.5000
39 -1147.6667 -4374.3333
40 -1392.5000 -1147.6667
41 4710.6667 -1392.5000
42 -2384.0000 4710.6667
43 2359.3333 -2384.0000
44 988.0000 2359.3333
45 -3767.0000 988.0000
46 2570.1667 -3767.0000
47 -1450.8000 2570.1667
48 5125.8333 -1450.8000
49 5465.5000 5125.8333
50 8647.6667 5465.5000
51 1858.3333 8647.6667
52 8007.5000 1858.3333
53 1766.6667 8007.5000
54 573.0000 1766.6667
55 2233.3333 573.0000
56 536.0000 2233.3333
57 867.0000 536.0000
58 1754.1667 867.0000
59 -1212.8000 1754.1667
60 2699.8333 -1212.8000
61 -258.5000 2699.8333
62 2830.6667 -258.5000
63 -1608.6667 2830.6667
64 4334.5000 -1608.6667
65 5622.6667 4334.5000
66 4928.0000 5622.6667
67 3129.3333 4928.0000
68 2203.0000 3129.3333
69 5555.0000 2203.0000
70 3466.1667 5555.0000
> 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/7614m1209843625.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/8wfqa1209843625.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/9bakw1209843625.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/103hhd1209843625.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/11z3v91209843625.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/12on441209843625.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/13pp2v1209843625.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/14zl5j1209843625.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/15pnod1209843625.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/16rtuz1209843625.tab")
+ }
>
> system("convert tmp/1hurs1209843625.ps tmp/1hurs1209843625.png")
> system("convert tmp/2xyeb1209843625.ps tmp/2xyeb1209843625.png")
> system("convert tmp/3zvk61209843625.ps tmp/3zvk61209843625.png")
> system("convert tmp/42ylf1209843625.ps tmp/42ylf1209843625.png")
> system("convert tmp/5la6f1209843625.ps tmp/5la6f1209843625.png")
> system("convert tmp/6y7yg1209843625.ps tmp/6y7yg1209843625.png")
> system("convert tmp/7614m1209843625.ps tmp/7614m1209843625.png")
> system("convert tmp/8wfqa1209843625.ps tmp/8wfqa1209843625.png")
> system("convert tmp/9bakw1209843625.ps tmp/9bakw1209843625.png")
> system("convert tmp/103hhd1209843625.ps tmp/103hhd1209843625.png")
>
>
> proc.time()
user system elapsed
2.830 1.605 3.733