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(8.1
+ ,10.9
+ ,115.6
+ ,92.9
+ ,7.7
+ ,10
+ ,127.1
+ ,107.7
+ ,7.5
+ ,9.2
+ ,123
+ ,103.5
+ ,7.6
+ ,9.2
+ ,122.2
+ ,91.1
+ ,7.8
+ ,9.5
+ ,126.4
+ ,79.8
+ ,7.8
+ ,9.6
+ ,112.7
+ ,71.9
+ ,7.8
+ ,9.5
+ ,105.8
+ ,82.9
+ ,7.5
+ ,9.1
+ ,120.9
+ ,90.1
+ ,7.5
+ ,8.9
+ ,116.3
+ ,100.7
+ ,7.1
+ ,9
+ ,115.7
+ ,90.7
+ ,7.5
+ ,10.1
+ ,127.9
+ ,108.8
+ ,7.5
+ ,10.3
+ ,108.3
+ ,44.1
+ ,7.6
+ ,10.2
+ ,121.1
+ ,93.6
+ ,7.7
+ ,9.6
+ ,128.6
+ ,107.4
+ ,7.7
+ ,9.2
+ ,123.1
+ ,96.5
+ ,7.9
+ ,9.3
+ ,127.7
+ ,93.6
+ ,8.1
+ ,9.4
+ ,126.6
+ ,76.5
+ ,8.2
+ ,9.4
+ ,118.4
+ ,76.7
+ ,8.2
+ ,9.2
+ ,110
+ ,84
+ ,8.2
+ ,9
+ ,129.6
+ ,103.3
+ ,7.9
+ ,9
+ ,115.8
+ ,88.5
+ ,7.3
+ ,9
+ ,125.9
+ ,99
+ ,6.9
+ ,9.8
+ ,128.4
+ ,105.9
+ ,6.6
+ ,10
+ ,114
+ ,44.7
+ ,6.7
+ ,9.8
+ ,125.6
+ ,94
+ ,6.9
+ ,9.3
+ ,128.5
+ ,107.1
+ ,7
+ ,9
+ ,136.6
+ ,104.8
+ ,7.1
+ ,9
+ ,133.1
+ ,102.5
+ ,7.2
+ ,9.1
+ ,124.6
+ ,77.7
+ ,7.1
+ ,9.1
+ ,123.5
+ ,85.2
+ ,6.9
+ ,9.1
+ ,117.2
+ ,91.3
+ ,7
+ ,9.2
+ ,135.5
+ ,106.5
+ ,6.8
+ ,8.8
+ ,124.8
+ ,92.4
+ ,6.4
+ ,8.3
+ ,127.8
+ ,97.5
+ ,6.7
+ ,8.4
+ ,133.1
+ ,107
+ ,6.6
+ ,8.1
+ ,125.7
+ ,51.1
+ ,6.4
+ ,7.7
+ ,128.4
+ ,98.6
+ ,6.3
+ ,7.9
+ ,131.9
+ ,102.2
+ ,6.2
+ ,7.9
+ ,146.3
+ ,114.3
+ ,6.5
+ ,8
+ ,140.6
+ ,99.4
+ ,6.8
+ ,7.9
+ ,129.5
+ ,72.5
+ ,6.8
+ ,7.6
+ ,132.4
+ ,92.3
+ ,6.4
+ ,7.1
+ ,125.9
+ ,99.4
+ ,6.1
+ ,6.8
+ ,126.9
+ ,85.9
+ ,5.8
+ ,6.5
+ ,135.8
+ ,109.4
+ ,6.1
+ ,6.9
+ ,129.5
+ ,97.6
+ ,7.2
+ ,8.2
+ ,130.2
+ ,104.7
+ ,7.3
+ ,8.7
+ ,133.8
+ ,56.9
+ ,6.9
+ ,8.3
+ ,123.3
+ ,86.7
+ ,6.1
+ ,7.9
+ ,140.7
+ ,108.5
+ ,5.8
+ ,7.5
+ ,145.9
+ ,103.4
+ ,6.2
+ ,7.8
+ ,128.5
+ ,86.2
+ ,7.1
+ ,8.3
+ ,135.9
+ ,71
+ ,7.7
+ ,8.4
+ ,120.2
+ ,75.9
+ ,7.9
+ ,8.2
+ ,119.2
+ ,87.1
+ ,7.7
+ ,7.7
+ ,132.5
+ ,102
+ ,7.4
+ ,7.2
+ ,130.5
+ ,88.5
+ ,7.5
+ ,7.3
+ ,124.8
+ ,87.8
+ ,8
+ ,8.1
+ ,136.7
+ ,100.8
+ ,8.1
+ ,8.5
+ ,129.2
+ ,50.6
+ ,8
+ ,8.4
+ ,127.9
+ ,85.9)
+ ,dim=c(4
+ ,61)
+ ,dimnames=list(c('mannen'
+ ,'vrouwen'
+ ,'voeding'
+ ,'bouw')
+ ,1:61))
> y <- array(NA,dim=c(4,61),dimnames=list(c('mannen','vrouwen','voeding','bouw'),1:61))
> 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 = '4'
> #'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
bouw mannen vrouwen voeding
1 92.9 8.1 10.9 115.6
2 107.7 7.7 10.0 127.1
3 103.5 7.5 9.2 123.0
4 91.1 7.6 9.2 122.2
5 79.8 7.8 9.5 126.4
6 71.9 7.8 9.6 112.7
7 82.9 7.8 9.5 105.8
8 90.1 7.5 9.1 120.9
9 100.7 7.5 8.9 116.3
10 90.7 7.1 9.0 115.7
11 108.8 7.5 10.1 127.9
12 44.1 7.5 10.3 108.3
13 93.6 7.6 10.2 121.1
14 107.4 7.7 9.6 128.6
15 96.5 7.7 9.2 123.1
16 93.6 7.9 9.3 127.7
17 76.5 8.1 9.4 126.6
18 76.7 8.2 9.4 118.4
19 84.0 8.2 9.2 110.0
20 103.3 8.2 9.0 129.6
21 88.5 7.9 9.0 115.8
22 99.0 7.3 9.0 125.9
23 105.9 6.9 9.8 128.4
24 44.7 6.6 10.0 114.0
25 94.0 6.7 9.8 125.6
26 107.1 6.9 9.3 128.5
27 104.8 7.0 9.0 136.6
28 102.5 7.1 9.0 133.1
29 77.7 7.2 9.1 124.6
30 85.2 7.1 9.1 123.5
31 91.3 6.9 9.1 117.2
32 106.5 7.0 9.2 135.5
33 92.4 6.8 8.8 124.8
34 97.5 6.4 8.3 127.8
35 107.0 6.7 8.4 133.1
36 51.1 6.6 8.1 125.7
37 98.6 6.4 7.7 128.4
38 102.2 6.3 7.9 131.9
39 114.3 6.2 7.9 146.3
40 99.4 6.5 8.0 140.6
41 72.5 6.8 7.9 129.5
42 92.3 6.8 7.6 132.4
43 99.4 6.4 7.1 125.9
44 85.9 6.1 6.8 126.9
45 109.4 5.8 6.5 135.8
46 97.6 6.1 6.9 129.5
47 104.7 7.2 8.2 130.2
48 56.9 7.3 8.7 133.8
49 86.7 6.9 8.3 123.3
50 108.5 6.1 7.9 140.7
51 103.4 5.8 7.5 145.9
52 86.2 6.2 7.8 128.5
53 71.0 7.1 8.3 135.9
54 75.9 7.7 8.4 120.2
55 87.1 7.9 8.2 119.2
56 102.0 7.7 7.7 132.5
57 88.5 7.4 7.2 130.5
58 87.8 7.5 7.3 124.8
59 100.8 8.0 8.1 136.7
60 50.6 8.1 8.5 129.2
61 85.9 8.0 8.4 127.9
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) mannen vrouwen voeding
-39.6100 -1.3797 2.1559 0.9589
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-40.823 -4.199 4.584 9.343 19.956
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -39.6100 54.8831 -0.722 0.47342
mannen -1.3797 3.6397 -0.379 0.70605
vrouwen 2.1559 2.6006 0.829 0.41057
voeding 0.9589 0.2792 3.434 0.00111 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 14.79 on 57 degrees of freedom
Multiple R-squared: 0.2158, Adjusted R-squared: 0.1746
F-statistic: 5.229 on 3 and 57 DF, p-value: 0.002938
> 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.08698089 0.1739618 0.9130191
[2,] 0.05583064 0.1116613 0.9441694
[3,] 0.05693880 0.1138776 0.9430612
[4,] 0.20036304 0.4007261 0.7996370
[5,] 0.12752403 0.2550481 0.8724760
[6,] 0.62011889 0.7597622 0.3798811
[7,] 0.52578118 0.9484376 0.4742188
[8,] 0.44527727 0.8905545 0.5547227
[9,] 0.36221351 0.7244270 0.6377865
[10,] 0.32229842 0.6445968 0.6777016
[11,] 0.37836421 0.7567284 0.6216358
[12,] 0.29612140 0.5922428 0.7038786
[13,] 0.29065787 0.5813157 0.7093421
[14,] 0.24507298 0.4901460 0.7549270
[15,] 0.20476981 0.4095396 0.7952302
[16,] 0.16433634 0.3286727 0.8356637
[17,] 0.13710829 0.2742166 0.8628917
[18,] 0.46460283 0.9292057 0.5353972
[19,] 0.38869622 0.7773924 0.6113038
[20,] 0.36710786 0.7342157 0.6328921
[21,] 0.32705960 0.6541192 0.6729404
[22,] 0.28622000 0.5724400 0.7137800
[23,] 0.27648878 0.5529776 0.7235112
[24,] 0.22125186 0.4425037 0.7787481
[25,] 0.20608018 0.4121604 0.7939198
[26,] 0.20071990 0.4014398 0.7992801
[27,] 0.17604637 0.3520927 0.8239536
[28,] 0.15304338 0.3060868 0.8469566
[29,] 0.17298835 0.3459767 0.8270116
[30,] 0.58959287 0.8208143 0.4104071
[31,] 0.53634861 0.9273028 0.4636514
[32,] 0.49965297 0.9993059 0.5003470
[33,] 0.47691294 0.9538259 0.5230871
[34,] 0.42670176 0.8534035 0.5732982
[35,] 0.47126499 0.9425300 0.5287350
[36,] 0.38863654 0.7772731 0.6113635
[37,] 0.33650499 0.6730100 0.6634950
[38,] 0.32636536 0.6527307 0.6736346
[39,] 0.26740723 0.5348145 0.7325928
[40,] 0.22278432 0.4455686 0.7772157
[41,] 0.29214313 0.5842863 0.7078569
[42,] 0.52355720 0.9528856 0.4764428
[43,] 0.43933878 0.8786776 0.5606612
[44,] 0.43299535 0.8659907 0.5670046
[45,] 0.34380809 0.6876162 0.6561919
[46,] 0.26973780 0.5394756 0.7302622
[47,] 0.21051321 0.4210264 0.7894868
[48,] 0.13133308 0.2626662 0.8686669
> postscript(file="/var/www/html/rcomp/tmp/1yka21258714448.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/2dpco1258714448.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/3f8271258714448.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/4sl6u1258714448.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/525r81258714448.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 = 61
Frequency = 1
1 2 3 4 5 6
9.3433048 14.5049450 15.6849953 4.1900431 -11.5079547 -6.4872952
7 8 9 10 11 12
11.3443584 4.5141681 19.9560539 9.7639057 14.3463420 -31.9913702
13 14 15 16 17 18
5.5889031 13.6290198 8.8650457 1.6146829 -14.3702341 -6.1696952
19 20 21 22 23 24
9.6158211 10.5535332 8.5717622 8.5595696 10.7858733 -37.4517627
25 26 27 28 29 30
1.2947184 12.9679259 3.6859699 4.8799131 -11.8474805 -3.4307130
31 32 33 34 35 36
8.4341076 6.0095300 2.7556413 5.5051577 10.1215675 -38.1741466
37 38 39 40 41 42
7.3233727 6.9982545 5.1528446 -4.0833940 -19.7106668 -2.0445696
43 44 45 46 47 48
11.8140231 -2.4119675 12.7871256 6.5794346 11.7232463 -40.4685841
49 50 51 52 53 54
-0.2901782 4.5844378 -5.0531360 -5.6640352 -27.7957547 -7.2295888
55 56 57 58 59 60
5.6363718 8.5856670 -2.3325982 2.3552278 2.9100494 -40.8229569
61
-4.1988318
> postscript(file="/var/www/html/rcomp/tmp/6i8tb1258714448.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 = 61
Frequency = 1
lag(myerror, k = 1) myerror
0 9.3433048 NA
1 14.5049450 9.3433048
2 15.6849953 14.5049450
3 4.1900431 15.6849953
4 -11.5079547 4.1900431
5 -6.4872952 -11.5079547
6 11.3443584 -6.4872952
7 4.5141681 11.3443584
8 19.9560539 4.5141681
9 9.7639057 19.9560539
10 14.3463420 9.7639057
11 -31.9913702 14.3463420
12 5.5889031 -31.9913702
13 13.6290198 5.5889031
14 8.8650457 13.6290198
15 1.6146829 8.8650457
16 -14.3702341 1.6146829
17 -6.1696952 -14.3702341
18 9.6158211 -6.1696952
19 10.5535332 9.6158211
20 8.5717622 10.5535332
21 8.5595696 8.5717622
22 10.7858733 8.5595696
23 -37.4517627 10.7858733
24 1.2947184 -37.4517627
25 12.9679259 1.2947184
26 3.6859699 12.9679259
27 4.8799131 3.6859699
28 -11.8474805 4.8799131
29 -3.4307130 -11.8474805
30 8.4341076 -3.4307130
31 6.0095300 8.4341076
32 2.7556413 6.0095300
33 5.5051577 2.7556413
34 10.1215675 5.5051577
35 -38.1741466 10.1215675
36 7.3233727 -38.1741466
37 6.9982545 7.3233727
38 5.1528446 6.9982545
39 -4.0833940 5.1528446
40 -19.7106668 -4.0833940
41 -2.0445696 -19.7106668
42 11.8140231 -2.0445696
43 -2.4119675 11.8140231
44 12.7871256 -2.4119675
45 6.5794346 12.7871256
46 11.7232463 6.5794346
47 -40.4685841 11.7232463
48 -0.2901782 -40.4685841
49 4.5844378 -0.2901782
50 -5.0531360 4.5844378
51 -5.6640352 -5.0531360
52 -27.7957547 -5.6640352
53 -7.2295888 -27.7957547
54 5.6363718 -7.2295888
55 8.5856670 5.6363718
56 -2.3325982 8.5856670
57 2.3552278 -2.3325982
58 2.9100494 2.3552278
59 -40.8229569 2.9100494
60 -4.1988318 -40.8229569
61 NA -4.1988318
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 14.5049450 9.3433048
[2,] 15.6849953 14.5049450
[3,] 4.1900431 15.6849953
[4,] -11.5079547 4.1900431
[5,] -6.4872952 -11.5079547
[6,] 11.3443584 -6.4872952
[7,] 4.5141681 11.3443584
[8,] 19.9560539 4.5141681
[9,] 9.7639057 19.9560539
[10,] 14.3463420 9.7639057
[11,] -31.9913702 14.3463420
[12,] 5.5889031 -31.9913702
[13,] 13.6290198 5.5889031
[14,] 8.8650457 13.6290198
[15,] 1.6146829 8.8650457
[16,] -14.3702341 1.6146829
[17,] -6.1696952 -14.3702341
[18,] 9.6158211 -6.1696952
[19,] 10.5535332 9.6158211
[20,] 8.5717622 10.5535332
[21,] 8.5595696 8.5717622
[22,] 10.7858733 8.5595696
[23,] -37.4517627 10.7858733
[24,] 1.2947184 -37.4517627
[25,] 12.9679259 1.2947184
[26,] 3.6859699 12.9679259
[27,] 4.8799131 3.6859699
[28,] -11.8474805 4.8799131
[29,] -3.4307130 -11.8474805
[30,] 8.4341076 -3.4307130
[31,] 6.0095300 8.4341076
[32,] 2.7556413 6.0095300
[33,] 5.5051577 2.7556413
[34,] 10.1215675 5.5051577
[35,] -38.1741466 10.1215675
[36,] 7.3233727 -38.1741466
[37,] 6.9982545 7.3233727
[38,] 5.1528446 6.9982545
[39,] -4.0833940 5.1528446
[40,] -19.7106668 -4.0833940
[41,] -2.0445696 -19.7106668
[42,] 11.8140231 -2.0445696
[43,] -2.4119675 11.8140231
[44,] 12.7871256 -2.4119675
[45,] 6.5794346 12.7871256
[46,] 11.7232463 6.5794346
[47,] -40.4685841 11.7232463
[48,] -0.2901782 -40.4685841
[49,] 4.5844378 -0.2901782
[50,] -5.0531360 4.5844378
[51,] -5.6640352 -5.0531360
[52,] -27.7957547 -5.6640352
[53,] -7.2295888 -27.7957547
[54,] 5.6363718 -7.2295888
[55,] 8.5856670 5.6363718
[56,] -2.3325982 8.5856670
[57,] 2.3552278 -2.3325982
[58,] 2.9100494 2.3552278
[59,] -40.8229569 2.9100494
[60,] -4.1988318 -40.8229569
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 14.5049450 9.3433048
2 15.6849953 14.5049450
3 4.1900431 15.6849953
4 -11.5079547 4.1900431
5 -6.4872952 -11.5079547
6 11.3443584 -6.4872952
7 4.5141681 11.3443584
8 19.9560539 4.5141681
9 9.7639057 19.9560539
10 14.3463420 9.7639057
11 -31.9913702 14.3463420
12 5.5889031 -31.9913702
13 13.6290198 5.5889031
14 8.8650457 13.6290198
15 1.6146829 8.8650457
16 -14.3702341 1.6146829
17 -6.1696952 -14.3702341
18 9.6158211 -6.1696952
19 10.5535332 9.6158211
20 8.5717622 10.5535332
21 8.5595696 8.5717622
22 10.7858733 8.5595696
23 -37.4517627 10.7858733
24 1.2947184 -37.4517627
25 12.9679259 1.2947184
26 3.6859699 12.9679259
27 4.8799131 3.6859699
28 -11.8474805 4.8799131
29 -3.4307130 -11.8474805
30 8.4341076 -3.4307130
31 6.0095300 8.4341076
32 2.7556413 6.0095300
33 5.5051577 2.7556413
34 10.1215675 5.5051577
35 -38.1741466 10.1215675
36 7.3233727 -38.1741466
37 6.9982545 7.3233727
38 5.1528446 6.9982545
39 -4.0833940 5.1528446
40 -19.7106668 -4.0833940
41 -2.0445696 -19.7106668
42 11.8140231 -2.0445696
43 -2.4119675 11.8140231
44 12.7871256 -2.4119675
45 6.5794346 12.7871256
46 11.7232463 6.5794346
47 -40.4685841 11.7232463
48 -0.2901782 -40.4685841
49 4.5844378 -0.2901782
50 -5.0531360 4.5844378
51 -5.6640352 -5.0531360
52 -27.7957547 -5.6640352
53 -7.2295888 -27.7957547
54 5.6363718 -7.2295888
55 8.5856670 5.6363718
56 -2.3325982 8.5856670
57 2.3552278 -2.3325982
58 2.9100494 2.3552278
59 -40.8229569 2.9100494
60 -4.1988318 -40.8229569
> 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/73rj21258714448.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/8td4h1258714448.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/9auas1258714448.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/10w2ca1258714448.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/1182e61258714448.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/121drc1258714448.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/13n7vj1258714448.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/14mavo1258714448.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/15x4cu1258714448.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/16d5ne1258714448.tab")
+ }
>
> system("convert tmp/1yka21258714448.ps tmp/1yka21258714448.png")
> system("convert tmp/2dpco1258714448.ps tmp/2dpco1258714448.png")
> system("convert tmp/3f8271258714448.ps tmp/3f8271258714448.png")
> system("convert tmp/4sl6u1258714448.ps tmp/4sl6u1258714448.png")
> system("convert tmp/525r81258714448.ps tmp/525r81258714448.png")
> system("convert tmp/6i8tb1258714448.ps tmp/6i8tb1258714448.png")
> system("convert tmp/73rj21258714448.ps tmp/73rj21258714448.png")
> system("convert tmp/8td4h1258714448.ps tmp/8td4h1258714448.png")
> system("convert tmp/9auas1258714448.ps tmp/9auas1258714448.png")
> system("convert tmp/10w2ca1258714448.ps tmp/10w2ca1258714448.png")
>
>
> proc.time()
user system elapsed
2.518 1.573 2.870