R version 2.9.0 (2009-04-17)
Copyright (C) 2009 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> x <- array(list(921365,0,987921,0,1132614,0,1332224,0,1418133,0,1411549,0,1695920,0,1636173,0,1539653,0,1395314,0,1127575,0,1036076,0,989236,0,1008380,0,1207763,0,1368839,0,1469798,0,1498721,0,1761769,0,1653214,0,1599104,0,1421179,0,1163995,0,1037735,0,1015407,0,1039210,0,1258049,0,1469445,0,1552346,0,1549144,0,1785895,0,1662335,0,1629440,0,1467430,0,1202209,0,1076982,0,1039367,1,1063449,1,1335135,1,1491602,1,1591972,1,1641248,1,1898849,1,1798580,1,1762444,1,1622044,1,1368955,1,1262973,1,1195650,1,1269530,1,1479279,1,1607819,1,1712466,1,1721766,1,1949843,1,1821326,1,1757802,1,1590367,1,1260647,1,1149235,1,1016367,1,1027885,1,1262159,1,1520854,1,1544144,1,1564709,1,1821776,1,1741365,1,1623386,1,1498658,1,1241822,1,1136029,1),dim=c(2,72),dimnames=list(c('bewegingen','dummy'),1:72))
> y <- array(NA,dim=c(2,72),dimnames=list(c('bewegingen','dummy'),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
bewegingen dummy M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 921365 0 1 0 0 0 0 0 0 0 0 0 0
2 987921 0 0 1 0 0 0 0 0 0 0 0 0
3 1132614 0 0 0 1 0 0 0 0 0 0 0 0
4 1332224 0 0 0 0 1 0 0 0 0 0 0 0
5 1418133 0 0 0 0 0 1 0 0 0 0 0 0
6 1411549 0 0 0 0 0 0 1 0 0 0 0 0
7 1695920 0 0 0 0 0 0 0 1 0 0 0 0
8 1636173 0 0 0 0 0 0 0 0 1 0 0 0
9 1539653 0 0 0 0 0 0 0 0 0 1 0 0
10 1395314 0 0 0 0 0 0 0 0 0 0 1 0
11 1127575 0 0 0 0 0 0 0 0 0 0 0 1
12 1036076 0 0 0 0 0 0 0 0 0 0 0 0
13 989236 0 1 0 0 0 0 0 0 0 0 0 0
14 1008380 0 0 1 0 0 0 0 0 0 0 0 0
15 1207763 0 0 0 1 0 0 0 0 0 0 0 0
16 1368839 0 0 0 0 1 0 0 0 0 0 0 0
17 1469798 0 0 0 0 0 1 0 0 0 0 0 0
18 1498721 0 0 0 0 0 0 1 0 0 0 0 0
19 1761769 0 0 0 0 0 0 0 1 0 0 0 0
20 1653214 0 0 0 0 0 0 0 0 1 0 0 0
21 1599104 0 0 0 0 0 0 0 0 0 1 0 0
22 1421179 0 0 0 0 0 0 0 0 0 0 1 0
23 1163995 0 0 0 0 0 0 0 0 0 0 0 1
24 1037735 0 0 0 0 0 0 0 0 0 0 0 0
25 1015407 0 1 0 0 0 0 0 0 0 0 0 0
26 1039210 0 0 1 0 0 0 0 0 0 0 0 0
27 1258049 0 0 0 1 0 0 0 0 0 0 0 0
28 1469445 0 0 0 0 1 0 0 0 0 0 0 0
29 1552346 0 0 0 0 0 1 0 0 0 0 0 0
30 1549144 0 0 0 0 0 0 1 0 0 0 0 0
31 1785895 0 0 0 0 0 0 0 1 0 0 0 0
32 1662335 0 0 0 0 0 0 0 0 1 0 0 0
33 1629440 0 0 0 0 0 0 0 0 0 1 0 0
34 1467430 0 0 0 0 0 0 0 0 0 0 1 0
35 1202209 0 0 0 0 0 0 0 0 0 0 0 1
36 1076982 0 0 0 0 0 0 0 0 0 0 0 0
37 1039367 1 1 0 0 0 0 0 0 0 0 0 0
38 1063449 1 0 1 0 0 0 0 0 0 0 0 0
39 1335135 1 0 0 1 0 0 0 0 0 0 0 0
40 1491602 1 0 0 0 1 0 0 0 0 0 0 0
41 1591972 1 0 0 0 0 1 0 0 0 0 0 0
42 1641248 1 0 0 0 0 0 1 0 0 0 0 0
43 1898849 1 0 0 0 0 0 0 1 0 0 0 0
44 1798580 1 0 0 0 0 0 0 0 1 0 0 0
45 1762444 1 0 0 0 0 0 0 0 0 1 0 0
46 1622044 1 0 0 0 0 0 0 0 0 0 1 0
47 1368955 1 0 0 0 0 0 0 0 0 0 0 1
48 1262973 1 0 0 0 0 0 0 0 0 0 0 0
49 1195650 1 1 0 0 0 0 0 0 0 0 0 0
50 1269530 1 0 1 0 0 0 0 0 0 0 0 0
51 1479279 1 0 0 1 0 0 0 0 0 0 0 0
52 1607819 1 0 0 0 1 0 0 0 0 0 0 0
53 1712466 1 0 0 0 0 1 0 0 0 0 0 0
54 1721766 1 0 0 0 0 0 1 0 0 0 0 0
55 1949843 1 0 0 0 0 0 0 1 0 0 0 0
56 1821326 1 0 0 0 0 0 0 0 1 0 0 0
57 1757802 1 0 0 0 0 0 0 0 0 1 0 0
58 1590367 1 0 0 0 0 0 0 0 0 0 1 0
59 1260647 1 0 0 0 0 0 0 0 0 0 0 1
60 1149235 1 0 0 0 0 0 0 0 0 0 0 0
61 1016367 1 1 0 0 0 0 0 0 0 0 0 0
62 1027885 1 0 1 0 0 0 0 0 0 0 0 0
63 1262159 1 0 0 1 0 0 0 0 0 0 0 0
64 1520854 1 0 0 0 1 0 0 0 0 0 0 0
65 1544144 1 0 0 0 0 1 0 0 0 0 0 0
66 1564709 1 0 0 0 0 0 1 0 0 0 0 0
67 1821776 1 0 0 0 0 0 0 1 0 0 0 0
68 1741365 1 0 0 0 0 0 0 0 1 0 0 0
69 1623386 1 0 0 0 0 0 0 0 0 1 0 0
70 1498658 1 0 0 0 0 0 0 0 0 0 1 0
71 1241822 1 0 0 0 0 0 0 0 0 0 0 1
72 1136029 1 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) dummy M1 M2 M3 M4
1048875 135260 -86940 -50442 162662 348626
M5 M6 M7 M8 M9 M10
431638 448018 702504 602327 535467 382660
M11
111029
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-105807.5 -44790.7 -972.7 41159.4 135837.5
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1048875 26330 39.835 < 2e-16 ***
dummy 135260 14606 9.261 4.21e-13 ***
M1 -86940 35776 -2.430 0.01816 *
M2 -50442 35776 -1.410 0.16380
M3 162662 35776 4.547 2.76e-05 ***
M4 348626 35776 9.745 6.71e-14 ***
M5 431638 35776 12.065 < 2e-16 ***
M6 448018 35776 12.523 < 2e-16 ***
M7 702504 35776 19.636 < 2e-16 ***
M8 602327 35776 16.836 < 2e-16 ***
M9 535466 35776 14.967 < 2e-16 ***
M10 382660 35776 10.696 1.95e-15 ***
M11 111029 35776 3.103 0.00294 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 61970 on 59 degrees of freedom
Multiple R-squared: 0.9565, Adjusted R-squared: 0.9476
F-statistic: 108 on 12 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.3080352483 0.6160704967 0.6919648
[2,] 0.2229708677 0.4459417354 0.7770291
[3,] 0.2470652299 0.4941304597 0.7529348
[4,] 0.2068998981 0.4137997962 0.7931001
[5,] 0.1271267484 0.2542534968 0.8728733
[6,] 0.0995672601 0.1991345202 0.9004327
[7,] 0.0616459788 0.1232919576 0.9383540
[8,] 0.0388708712 0.0777417424 0.9611291
[9,] 0.0211445598 0.0422891195 0.9788554
[10,] 0.0189675437 0.0379350873 0.9810325
[11,] 0.0125073633 0.0250147266 0.9874926
[12,] 0.0184042265 0.0368084529 0.9815958
[13,] 0.0424809905 0.0849619810 0.9575190
[14,] 0.0620406402 0.1240812804 0.9379594
[15,] 0.0671784588 0.1343569177 0.9328215
[16,] 0.0514540320 0.1029080641 0.9485460
[17,] 0.0328674284 0.0657348568 0.9671326
[18,] 0.0249575116 0.0499150232 0.9750425
[19,] 0.0184724689 0.0369449378 0.9815275
[20,] 0.0132029457 0.0264058913 0.9867971
[21,] 0.0082794508 0.0165589016 0.9917205
[22,] 0.0050567728 0.0101135456 0.9949432
[23,] 0.0032992605 0.0065985209 0.9967007
[24,] 0.0025190053 0.0050380105 0.9974810
[25,] 0.0015364459 0.0030728918 0.9984636
[26,] 0.0008181443 0.0016362887 0.9991819
[27,] 0.0005236032 0.0010472064 0.9994764
[28,] 0.0002912049 0.0005824099 0.9997088
[29,] 0.0001470336 0.0002940673 0.9998530
[30,] 0.0001045158 0.0002090316 0.9998955
[31,] 0.0000897854 0.0001795708 0.9999102
[32,] 0.0001119616 0.0002239232 0.9998880
[33,] 0.0001476566 0.0002953132 0.9998523
[34,] 0.0004412579 0.0008825159 0.9995587
[35,] 0.0059668299 0.0119336597 0.9940332
[36,] 0.0410767539 0.0821535079 0.9589232
[37,] 0.0364467794 0.0728935587 0.9635532
[38,] 0.0850917660 0.1701835321 0.9149082
[39,] 0.1771150211 0.3542300421 0.8228850
[40,] 0.2573431876 0.5146863752 0.7426568
[41,] 0.2222005720 0.4444011440 0.7777994
> postscript(file="/var/www/html/rcomp/tmp/1m1sk1293456591.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/www/html/rcomp/tmp/2m1sk1293456591.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/www/html/rcomp/tmp/3m1sk1293456591.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/www/html/rcomp/tmp/4xsrn1293456591.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/www/html/rcomp/tmp/5xsrn1293456591.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 = 72
Frequency = 1
1 2 3 4 5 6
-40570.333 -10511.500 -78922.500 -65276.500 -62380.167 -85343.833
7 8 9 10 11 12
-55458.667 -15029.167 -44688.500 -36221.333 -32328.833 -12799.000
13 14 15 16 17 18
27300.667 9947.500 -3773.500 -28661.500 -10715.167 1828.167
19 20 21 22 23 24
10390.333 2011.833 14762.500 -10356.333 4091.167 -11140.000
25 26 27 28 29 30
53471.667 40777.500 46512.500 71944.500 71832.833 52251.167
31 32 33 34 35 36
34516.333 11132.833 45098.500 35894.667 42305.167 28107.000
37 38 39 40 41 42
-57828.333 -70243.500 -11661.500 -41158.500 -23801.167 9095.167
43 44 45 46 47 48
12210.333 12117.833 42842.500 55248.667 73791.167 78838.000
49 50 51 52 53 54
98454.667 135837.500 132482.500 75058.500 96692.833 89613.167
55 56 57 58 59 60
63204.333 34863.833 38200.500 23571.667 -34516.833 -34900.000
61 62 63 64 65 66
-80828.333 -105807.500 -84637.500 -11906.500 -71629.167 -67443.833
67 68 69 70 71 72
-64862.667 -45097.167 -96215.500 -68137.333 -53341.833 -48106.000
> postscript(file="/var/www/html/rcomp/tmp/6xsrn1293456591.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 = 72
Frequency = 1
lag(myerror, k = 1) myerror
0 -40570.333 NA
1 -10511.500 -40570.333
2 -78922.500 -10511.500
3 -65276.500 -78922.500
4 -62380.167 -65276.500
5 -85343.833 -62380.167
6 -55458.667 -85343.833
7 -15029.167 -55458.667
8 -44688.500 -15029.167
9 -36221.333 -44688.500
10 -32328.833 -36221.333
11 -12799.000 -32328.833
12 27300.667 -12799.000
13 9947.500 27300.667
14 -3773.500 9947.500
15 -28661.500 -3773.500
16 -10715.167 -28661.500
17 1828.167 -10715.167
18 10390.333 1828.167
19 2011.833 10390.333
20 14762.500 2011.833
21 -10356.333 14762.500
22 4091.167 -10356.333
23 -11140.000 4091.167
24 53471.667 -11140.000
25 40777.500 53471.667
26 46512.500 40777.500
27 71944.500 46512.500
28 71832.833 71944.500
29 52251.167 71832.833
30 34516.333 52251.167
31 11132.833 34516.333
32 45098.500 11132.833
33 35894.667 45098.500
34 42305.167 35894.667
35 28107.000 42305.167
36 -57828.333 28107.000
37 -70243.500 -57828.333
38 -11661.500 -70243.500
39 -41158.500 -11661.500
40 -23801.167 -41158.500
41 9095.167 -23801.167
42 12210.333 9095.167
43 12117.833 12210.333
44 42842.500 12117.833
45 55248.667 42842.500
46 73791.167 55248.667
47 78838.000 73791.167
48 98454.667 78838.000
49 135837.500 98454.667
50 132482.500 135837.500
51 75058.500 132482.500
52 96692.833 75058.500
53 89613.167 96692.833
54 63204.333 89613.167
55 34863.833 63204.333
56 38200.500 34863.833
57 23571.667 38200.500
58 -34516.833 23571.667
59 -34900.000 -34516.833
60 -80828.333 -34900.000
61 -105807.500 -80828.333
62 -84637.500 -105807.500
63 -11906.500 -84637.500
64 -71629.167 -11906.500
65 -67443.833 -71629.167
66 -64862.667 -67443.833
67 -45097.167 -64862.667
68 -96215.500 -45097.167
69 -68137.333 -96215.500
70 -53341.833 -68137.333
71 -48106.000 -53341.833
72 NA -48106.000
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -10511.500 -40570.333
[2,] -78922.500 -10511.500
[3,] -65276.500 -78922.500
[4,] -62380.167 -65276.500
[5,] -85343.833 -62380.167
[6,] -55458.667 -85343.833
[7,] -15029.167 -55458.667
[8,] -44688.500 -15029.167
[9,] -36221.333 -44688.500
[10,] -32328.833 -36221.333
[11,] -12799.000 -32328.833
[12,] 27300.667 -12799.000
[13,] 9947.500 27300.667
[14,] -3773.500 9947.500
[15,] -28661.500 -3773.500
[16,] -10715.167 -28661.500
[17,] 1828.167 -10715.167
[18,] 10390.333 1828.167
[19,] 2011.833 10390.333
[20,] 14762.500 2011.833
[21,] -10356.333 14762.500
[22,] 4091.167 -10356.333
[23,] -11140.000 4091.167
[24,] 53471.667 -11140.000
[25,] 40777.500 53471.667
[26,] 46512.500 40777.500
[27,] 71944.500 46512.500
[28,] 71832.833 71944.500
[29,] 52251.167 71832.833
[30,] 34516.333 52251.167
[31,] 11132.833 34516.333
[32,] 45098.500 11132.833
[33,] 35894.667 45098.500
[34,] 42305.167 35894.667
[35,] 28107.000 42305.167
[36,] -57828.333 28107.000
[37,] -70243.500 -57828.333
[38,] -11661.500 -70243.500
[39,] -41158.500 -11661.500
[40,] -23801.167 -41158.500
[41,] 9095.167 -23801.167
[42,] 12210.333 9095.167
[43,] 12117.833 12210.333
[44,] 42842.500 12117.833
[45,] 55248.667 42842.500
[46,] 73791.167 55248.667
[47,] 78838.000 73791.167
[48,] 98454.667 78838.000
[49,] 135837.500 98454.667
[50,] 132482.500 135837.500
[51,] 75058.500 132482.500
[52,] 96692.833 75058.500
[53,] 89613.167 96692.833
[54,] 63204.333 89613.167
[55,] 34863.833 63204.333
[56,] 38200.500 34863.833
[57,] 23571.667 38200.500
[58,] -34516.833 23571.667
[59,] -34900.000 -34516.833
[60,] -80828.333 -34900.000
[61,] -105807.500 -80828.333
[62,] -84637.500 -105807.500
[63,] -11906.500 -84637.500
[64,] -71629.167 -11906.500
[65,] -67443.833 -71629.167
[66,] -64862.667 -67443.833
[67,] -45097.167 -64862.667
[68,] -96215.500 -45097.167
[69,] -68137.333 -96215.500
[70,] -53341.833 -68137.333
[71,] -48106.000 -53341.833
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -10511.500 -40570.333
2 -78922.500 -10511.500
3 -65276.500 -78922.500
4 -62380.167 -65276.500
5 -85343.833 -62380.167
6 -55458.667 -85343.833
7 -15029.167 -55458.667
8 -44688.500 -15029.167
9 -36221.333 -44688.500
10 -32328.833 -36221.333
11 -12799.000 -32328.833
12 27300.667 -12799.000
13 9947.500 27300.667
14 -3773.500 9947.500
15 -28661.500 -3773.500
16 -10715.167 -28661.500
17 1828.167 -10715.167
18 10390.333 1828.167
19 2011.833 10390.333
20 14762.500 2011.833
21 -10356.333 14762.500
22 4091.167 -10356.333
23 -11140.000 4091.167
24 53471.667 -11140.000
25 40777.500 53471.667
26 46512.500 40777.500
27 71944.500 46512.500
28 71832.833 71944.500
29 52251.167 71832.833
30 34516.333 52251.167
31 11132.833 34516.333
32 45098.500 11132.833
33 35894.667 45098.500
34 42305.167 35894.667
35 28107.000 42305.167
36 -57828.333 28107.000
37 -70243.500 -57828.333
38 -11661.500 -70243.500
39 -41158.500 -11661.500
40 -23801.167 -41158.500
41 9095.167 -23801.167
42 12210.333 9095.167
43 12117.833 12210.333
44 42842.500 12117.833
45 55248.667 42842.500
46 73791.167 55248.667
47 78838.000 73791.167
48 98454.667 78838.000
49 135837.500 98454.667
50 132482.500 135837.500
51 75058.500 132482.500
52 96692.833 75058.500
53 89613.167 96692.833
54 63204.333 89613.167
55 34863.833 63204.333
56 38200.500 34863.833
57 23571.667 38200.500
58 -34516.833 23571.667
59 -34900.000 -34516.833
60 -80828.333 -34900.000
61 -105807.500 -80828.333
62 -84637.500 -105807.500
63 -11906.500 -84637.500
64 -71629.167 -11906.500
65 -67443.833 -71629.167
66 -64862.667 -67443.833
67 -45097.167 -64862.667
68 -96215.500 -45097.167
69 -68137.333 -96215.500
70 -53341.833 -68137.333
71 -48106.000 -53341.833
> 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/7pjq81293456591.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/www/html/rcomp/tmp/8ibqt1293456591.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/www/html/rcomp/tmp/9ibqt1293456591.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')
hat values (leverages) are all = 0.1805556
and there are no factor predictors; no plot no. 5
> par(opar)
> dev.off()
null device
1
> if (n > n25) {
+ postscript(file="/var/www/html/rcomp/tmp/10ibqt1293456591.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/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/117v8f1293456592.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/12hm7i1293456592.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/136nmb1293456592.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/1495lh1293456592.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/152xkk1293456592.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/16y70b1293456592.tab")
+ }
>
> try(system("convert tmp/1m1sk1293456591.ps tmp/1m1sk1293456591.png",intern=TRUE))
character(0)
> try(system("convert tmp/2m1sk1293456591.ps tmp/2m1sk1293456591.png",intern=TRUE))
character(0)
> try(system("convert tmp/3m1sk1293456591.ps tmp/3m1sk1293456591.png",intern=TRUE))
character(0)
> try(system("convert tmp/4xsrn1293456591.ps tmp/4xsrn1293456591.png",intern=TRUE))
character(0)
> try(system("convert tmp/5xsrn1293456591.ps tmp/5xsrn1293456591.png",intern=TRUE))
character(0)
> try(system("convert tmp/6xsrn1293456591.ps tmp/6xsrn1293456591.png",intern=TRUE))
character(0)
> try(system("convert tmp/7pjq81293456591.ps tmp/7pjq81293456591.png",intern=TRUE))
character(0)
> try(system("convert tmp/8ibqt1293456591.ps tmp/8ibqt1293456591.png",intern=TRUE))
character(0)
> try(system("convert tmp/9ibqt1293456591.ps tmp/9ibqt1293456591.png",intern=TRUE))
character(0)
> try(system("convert tmp/10ibqt1293456591.ps tmp/10ibqt1293456591.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
2.587 1.642 6.051