R version 2.13.0 (2011-04-13)
Copyright (C) 2011 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: i486-pc-linux-gnu (32-bit)
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(6654000
+ ,5712000
+ ,3.3
+ ,38.6
+ ,645
+ ,3
+ ,5
+ ,3
+ ,1000
+ ,6600
+ ,8.3
+ ,4.5
+ ,42
+ ,3
+ ,1
+ ,3
+ ,3385
+ ,44500
+ ,12.5
+ ,14
+ ,60
+ ,1
+ ,1
+ ,1
+ ,0.92
+ ,5700
+ ,16.5
+ ,0
+ ,25
+ ,5
+ ,2
+ ,3
+ ,2547000
+ ,4603000
+ ,3.9
+ ,69
+ ,624
+ ,3
+ ,5
+ ,4
+ ,10550
+ ,179500
+ ,9.8
+ ,27
+ ,180
+ ,4
+ ,4
+ ,4
+ ,0.023
+ ,0.3
+ ,19.7
+ ,19
+ ,35
+ ,1
+ ,1
+ ,1
+ ,160000
+ ,169000
+ ,6.2
+ ,30.4
+ ,392
+ ,4
+ ,5
+ ,4
+ ,3300
+ ,25600
+ ,14.5
+ ,28
+ ,63
+ ,1
+ ,2
+ ,1
+ ,52160
+ ,440000
+ ,9.7
+ ,50
+ ,230
+ ,1
+ ,1
+ ,1
+ ,0.425
+ ,6400
+ ,12.5
+ ,7
+ ,112
+ ,5
+ ,4
+ ,4
+ ,465000
+ ,423000
+ ,3.9
+ ,30
+ ,281
+ ,5
+ ,5
+ ,5
+ ,0.55
+ ,2400
+ ,10.3
+ ,0
+ ,0
+ ,2
+ ,1
+ ,2
+ ,187100
+ ,419000
+ ,3.1
+ ,40
+ ,365
+ ,5
+ ,5
+ ,5
+ ,0.075
+ ,1200
+ ,8.4
+ ,3.5
+ ,42
+ ,1
+ ,1
+ ,1
+ ,3000
+ ,25000
+ ,8.6
+ ,50
+ ,28
+ ,2
+ ,2
+ ,2
+ ,0.785
+ ,3500
+ ,10.7
+ ,6
+ ,42
+ ,2
+ ,2
+ ,2
+ ,0.2
+ ,5000
+ ,10.7
+ ,10.4
+ ,120
+ ,2
+ ,2
+ ,2
+ ,1410
+ ,17500
+ ,6.1
+ ,34
+ ,0
+ ,1
+ ,2
+ ,1
+ ,60000
+ ,81000
+ ,18.1
+ ,7
+ ,0
+ ,1
+ ,1
+ ,1
+ ,27660
+ ,115000
+ ,3.8
+ ,20
+ ,148
+ ,5
+ ,5
+ ,5
+ ,0.12
+ ,1000
+ ,14.4
+ ,3.9
+ ,16
+ ,3
+ ,1
+ ,2
+ ,207000
+ ,406000
+ ,12
+ ,39.3
+ ,252
+ ,1
+ ,4
+ ,1
+ ,85000
+ ,325000
+ ,6.2
+ ,41
+ ,310
+ ,1
+ ,3
+ ,1
+ ,36330
+ ,119500
+ ,13
+ ,16.2
+ ,63
+ ,1
+ ,1
+ ,1
+ ,0.101
+ ,4000
+ ,13.8
+ ,9
+ ,28
+ ,5
+ ,1
+ ,3
+ ,1040
+ ,5500
+ ,8.2
+ ,7.6
+ ,68
+ ,5
+ ,3
+ ,4
+ ,521000
+ ,655000
+ ,2.9
+ ,46
+ ,336
+ ,5
+ ,5
+ ,5
+ ,100000
+ ,157000
+ ,10.8
+ ,22.4
+ ,100
+ ,1
+ ,1
+ ,1
+ ,0.005
+ ,0.14
+ ,9.1
+ ,2.6
+ ,21.5
+ ,5
+ ,2
+ ,4
+ ,0.01
+ ,0.25
+ ,19.9
+ ,24
+ ,50
+ ,1
+ ,1
+ ,1
+ ,62000
+ ,1320000
+ ,8
+ ,100
+ ,267
+ ,1
+ ,1
+ ,1
+ ,0.122
+ ,3000
+ ,10.6
+ ,0
+ ,30
+ ,2
+ ,1
+ ,1
+ ,1350
+ ,8100
+ ,11.2
+ ,0
+ ,45
+ ,3
+ ,1
+ ,3
+ ,0.023
+ ,0.4
+ ,13.2
+ ,3.2
+ ,19
+ ,4
+ ,1
+ ,3
+ ,0.048
+ ,0.33
+ ,12.8
+ ,2
+ ,30
+ ,4
+ ,1
+ ,3
+ ,1700
+ ,6300
+ ,19.4
+ ,5
+ ,12
+ ,2
+ ,1
+ ,1
+ ,3500
+ ,10800
+ ,17.4
+ ,6.5
+ ,120
+ ,2
+ ,1
+ ,1
+ ,0.48
+ ,15500
+ ,17
+ ,12
+ ,140
+ ,2
+ ,2
+ ,2
+ ,10000
+ ,115000
+ ,10.9
+ ,20.2
+ ,170
+ ,4
+ ,4
+ ,4
+ ,1620
+ ,11400
+ ,13.7
+ ,13
+ ,17
+ ,2
+ ,1
+ ,2
+ ,192000
+ ,180000
+ ,8.4
+ ,27
+ ,115
+ ,4
+ ,4
+ ,4
+ ,2500
+ ,12100
+ ,8.4
+ ,18
+ ,31
+ ,5
+ ,5
+ ,5
+ ,4288
+ ,39200
+ ,12.5
+ ,13.7
+ ,63
+ ,2
+ ,2
+ ,2
+ ,0.28
+ ,1900
+ ,13.2
+ ,4.7
+ ,21
+ ,3
+ ,1
+ ,3
+ ,4235
+ ,50400
+ ,9.8
+ ,9.8
+ ,52
+ ,1
+ ,1
+ ,1
+ ,6800
+ ,179000
+ ,9.6
+ ,29
+ ,164
+ ,2
+ ,3
+ ,2
+ ,0.75
+ ,12300
+ ,6.6
+ ,7
+ ,225
+ ,2
+ ,2
+ ,2
+ ,3600
+ ,21000
+ ,5.4
+ ,6
+ ,225
+ ,3
+ ,2
+ ,3
+ ,14830
+ ,98200
+ ,2.6
+ ,17
+ ,150
+ ,5
+ ,5
+ ,5
+ ,55500
+ ,175000
+ ,3.8
+ ,20
+ ,151
+ ,5
+ ,5
+ ,5
+ ,1400
+ ,12500
+ ,11
+ ,12.7
+ ,90
+ ,2
+ ,2
+ ,2
+ ,0.06
+ ,1000
+ ,10.3
+ ,3.5
+ ,0
+ ,3
+ ,1
+ ,2
+ ,0.9
+ ,2600
+ ,13.3
+ ,4.5
+ ,60
+ ,2
+ ,1
+ ,2
+ ,2000
+ ,12300
+ ,5.4
+ ,7.5
+ ,200
+ ,3
+ ,1
+ ,3
+ ,0.104
+ ,2500
+ ,15.8
+ ,2.3
+ ,46
+ ,3
+ ,2
+ ,2
+ ,4190
+ ,58000
+ ,10.3
+ ,24
+ ,210
+ ,4
+ ,3
+ ,4
+ ,3500
+ ,3900
+ ,19.4
+ ,3
+ ,14
+ ,2
+ ,1
+ ,1)
+ ,dim=c(8
+ ,58)
+ ,dimnames=list(c('gewicht'
+ ,'brein'
+ ,'totaleslaap'
+ ,'levensduur'
+ ,'drachttijd'
+ ,'jager?'
+ ,'blootgesteldheidslaap'
+ ,'algemeengevaar')
+ ,1:58))
> y <- array(NA,dim=c(8,58),dimnames=list(c('gewicht','brein','totaleslaap','levensduur','drachttijd','jager?','blootgesteldheidslaap','algemeengevaar'),1:58))
> 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 = '3'
> #'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
totaleslaap gewicht brein levensduur drachttijd jager?
1 3.3 6.654e+06 5.712e+06 38.6 645.0 3
2 8.3 1.000e+03 6.600e+03 4.5 42.0 3
3 12.5 3.385e+03 4.450e+04 14.0 60.0 1
4 16.5 9.200e-01 5.700e+03 0.0 25.0 5
5 3.9 2.547e+06 4.603e+06 69.0 624.0 3
6 9.8 1.055e+04 1.795e+05 27.0 180.0 4
7 19.7 2.300e-02 3.000e-01 19.0 35.0 1
8 6.2 1.600e+05 1.690e+05 30.4 392.0 4
9 14.5 3.300e+03 2.560e+04 28.0 63.0 1
10 9.7 5.216e+04 4.400e+05 50.0 230.0 1
11 12.5 4.250e-01 6.400e+03 7.0 112.0 5
12 3.9 4.650e+05 4.230e+05 30.0 281.0 5
13 10.3 5.500e-01 2.400e+03 0.0 0.0 2
14 3.1 1.871e+05 4.190e+05 40.0 365.0 5
15 8.4 7.500e-02 1.200e+03 3.5 42.0 1
16 8.6 3.000e+03 2.500e+04 50.0 28.0 2
17 10.7 7.850e-01 3.500e+03 6.0 42.0 2
18 10.7 2.000e-01 5.000e+03 10.4 120.0 2
19 6.1 1.410e+03 1.750e+04 34.0 0.0 1
20 18.1 6.000e+04 8.100e+04 7.0 0.0 1
21 3.8 2.766e+04 1.150e+05 20.0 148.0 5
22 14.4 1.200e-01 1.000e+03 3.9 16.0 3
23 12.0 2.070e+05 4.060e+05 39.3 252.0 1
24 6.2 8.500e+04 3.250e+05 41.0 310.0 1
25 13.0 3.633e+04 1.195e+05 16.2 63.0 1
26 13.8 1.010e-01 4.000e+03 9.0 28.0 5
27 8.2 1.040e+03 5.500e+03 7.6 68.0 5
28 2.9 5.210e+05 6.550e+05 46.0 336.0 5
29 10.8 1.000e+05 1.570e+05 22.4 100.0 1
30 9.1 5.000e-03 1.400e-01 2.6 21.5 5
31 19.9 1.000e-02 2.500e-01 24.0 50.0 1
32 8.0 6.200e+04 1.320e+06 100.0 267.0 1
33 10.6 1.220e-01 3.000e+03 0.0 30.0 2
34 11.2 1.350e+03 8.100e+03 0.0 45.0 3
35 13.2 2.300e-02 4.000e-01 3.2 19.0 4
36 12.8 4.800e-02 3.300e-01 2.0 30.0 4
37 19.4 1.700e+03 6.300e+03 5.0 12.0 2
38 17.4 3.500e+03 1.080e+04 6.5 120.0 2
39 17.0 4.800e-01 1.550e+04 12.0 140.0 2
40 10.9 1.000e+04 1.150e+05 20.2 170.0 4
41 13.7 1.620e+03 1.140e+04 13.0 17.0 2
42 8.4 1.920e+05 1.800e+05 27.0 115.0 4
43 8.4 2.500e+03 1.210e+04 18.0 31.0 5
44 12.5 4.288e+03 3.920e+04 13.7 63.0 2
45 13.2 2.800e-01 1.900e+03 4.7 21.0 3
46 9.8 4.235e+03 5.040e+04 9.8 52.0 1
47 9.6 6.800e+03 1.790e+05 29.0 164.0 2
48 6.6 7.500e-01 1.230e+04 7.0 225.0 2
49 5.4 3.600e+03 2.100e+04 6.0 225.0 3
50 2.6 1.483e+04 9.820e+04 17.0 150.0 5
51 3.8 5.550e+04 1.750e+05 20.0 151.0 5
52 11.0 1.400e+03 1.250e+04 12.7 90.0 2
53 10.3 6.000e-02 1.000e+03 3.5 0.0 3
54 13.3 9.000e-01 2.600e+03 4.5 60.0 2
55 5.4 2.000e+03 1.230e+04 7.5 200.0 3
56 15.8 1.040e-01 2.500e+03 2.3 46.0 3
57 10.3 4.190e+03 5.800e+04 24.0 210.0 4
58 19.4 3.500e+03 3.900e+03 3.0 14.0 2
blootgesteldheidslaap algemeengevaar
1 5 3
2 1 3
3 1 1
4 2 3
5 5 4
6 4 4
7 1 1
8 5 4
9 2 1
10 1 1
11 4 4
12 5 5
13 1 2
14 5 5
15 1 1
16 2 2
17 2 2
18 2 2
19 2 1
20 1 1
21 5 5
22 1 2
23 4 1
24 3 1
25 1 1
26 1 3
27 3 4
28 5 5
29 1 1
30 2 4
31 1 1
32 1 1
33 1 1
34 1 3
35 1 3
36 1 3
37 1 1
38 1 1
39 2 2
40 4 4
41 1 2
42 4 4
43 5 5
44 2 2
45 1 3
46 1 1
47 3 2
48 2 2
49 2 3
50 5 5
51 5 5
52 2 2
53 1 2
54 1 2
55 1 3
56 2 2
57 3 4
58 1 1
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) gewicht brein
1.590e+01 -1.612e-06 2.073e-06
levensduur drachttijd `jager?`
-4.771e-02 -1.199e-02 2.075e+00
blootgesteldheidslaap algemeengevaar
3.131e-01 -3.848e+00
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-7.0693 -1.9958 -0.0133 1.7292 7.2004
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.590e+01 1.068e+00 14.885 < 2e-16 ***
gewicht -1.612e-06 1.571e-06 -1.026 0.30999
brein 2.073e-06 1.727e-06 1.200 0.23565
levensduur -4.771e-02 3.713e-02 -1.285 0.20473
drachttijd -1.199e-02 6.433e-03 -1.864 0.06815 .
`jager?` 2.075e+00 8.919e-01 2.326 0.02409 *
blootgesteldheidslaap 3.131e-01 5.801e-01 0.540 0.59179
algemeengevaar -3.848e+00 1.108e+00 -3.475 0.00107 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.114 on 50 degrees of freedom
Multiple R-squared: 0.5993, Adjusted R-squared: 0.5432
F-statistic: 10.68 on 7 and 50 DF, p-value: 3.88e-08
> 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.6147559 0.7704882 0.3852441
[2,] 0.4674998 0.9349997 0.5325002
[3,] 0.4062865 0.8125730 0.5937135
[4,] 0.2781160 0.5562321 0.7218840
[5,] 0.4495186 0.8990371 0.5504814
[6,] 0.5349282 0.9301436 0.4650718
[7,] 0.4439301 0.8878602 0.5560699
[8,] 0.3392574 0.6785148 0.6607426
[9,] 0.8183199 0.3633602 0.1816801
[10,] 0.8727106 0.2545788 0.1272894
[11,] 0.8305161 0.3389678 0.1694839
[12,] 0.7658252 0.4683497 0.2341748
[13,] 0.7050570 0.5898860 0.2949430
[14,] 0.7772059 0.4455882 0.2227941
[15,] 0.7095072 0.5809856 0.2904928
[16,] 0.6591560 0.6816880 0.3408440
[17,] 0.6346220 0.7307560 0.3653780
[18,] 0.6043999 0.7912002 0.3956001
[19,] 0.5617940 0.8764119 0.4382060
[20,] 0.5057179 0.9885642 0.4942821
[21,] 0.6624095 0.6751809 0.3375905
[22,] 0.5910290 0.8179420 0.4089710
[23,] 0.7188639 0.5622723 0.2811361
[24,] 0.6858619 0.6282761 0.3141381
[25,] 0.6023932 0.7952135 0.3976068
[26,] 0.5104331 0.9791339 0.4895669
[27,] 0.4741982 0.9483964 0.5258018
[28,] 0.4063456 0.8126912 0.5936544
[29,] 0.6313684 0.7372632 0.3686316
[30,] 0.7311427 0.5377146 0.2688573
[31,] 0.6401059 0.7197882 0.3598941
[32,] 0.5454031 0.9091937 0.4545969
[33,] 0.4529240 0.9058480 0.5470760
[34,] 0.3337699 0.6675397 0.6662301
[35,] 0.3524589 0.7049179 0.6475411
[36,] 0.3026858 0.6053716 0.6973142
[37,] 0.1939764 0.3879529 0.8060236
> postscript(file="/var/wessaorg/rcomp/tmp/1gm8e1321988471.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/wessaorg/rcomp/tmp/227yk1321988471.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/wessaorg/rcomp/tmp/32q291321988471.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/wessaorg/rcomp/tmp/4sqjk1321988471.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/wessaorg/rcomp/tmp/5s7ak1321988471.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 = 58
Frequency = 1
1 2 3 4 5 6
-0.389998991 -1.891414074 -0.643598441 1.427259050 0.936508116 2.828172527
7 8 9 10 11 12
6.581913218 1.882652063 1.786283920 -0.428332015 2.025278451 -0.029471374
13 14 15 16 17 18
-2.376073504 0.215528546 -5.376151380 -1.709777753 -1.501448029 -0.359094720
19 20 21 22 23 24
-7.069319050 3.918355672 -2.268128703 0.029875362 1.005730484 -3.733097998
25 26 27 28 29 30
-0.105041858 -0.490714341 -2.457174786 0.002862262 -1.540574800 -2.030615255
31 32 33 34 35 36
7.200383000 -1.107525711 -5.565736347 0.827318374 0.607835731 0.282515485
37 38 39 40 41 42
3.252829771 2.613320646 6.235355970 3.616627163 1.832031098 0.939980135
43 44 45 46 47 48
1.005932161 0.850704106 2.774387270 -3.650802594 -0.706810160 -3.377079716
49 50 51 52 53 54
-2.863722988 -3.573125462 -2.311668241 -0.322469968 -4.281113406 1.557855816
55 56 57 58
-2.763432878 1.397132327 4.099615099 3.189270718
> postscript(file="/var/wessaorg/rcomp/tmp/6fet31321988471.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 = 58
Frequency = 1
lag(myerror, k = 1) myerror
0 -0.389998991 NA
1 -1.891414074 -0.389998991
2 -0.643598441 -1.891414074
3 1.427259050 -0.643598441
4 0.936508116 1.427259050
5 2.828172527 0.936508116
6 6.581913218 2.828172527
7 1.882652063 6.581913218
8 1.786283920 1.882652063
9 -0.428332015 1.786283920
10 2.025278451 -0.428332015
11 -0.029471374 2.025278451
12 -2.376073504 -0.029471374
13 0.215528546 -2.376073504
14 -5.376151380 0.215528546
15 -1.709777753 -5.376151380
16 -1.501448029 -1.709777753
17 -0.359094720 -1.501448029
18 -7.069319050 -0.359094720
19 3.918355672 -7.069319050
20 -2.268128703 3.918355672
21 0.029875362 -2.268128703
22 1.005730484 0.029875362
23 -3.733097998 1.005730484
24 -0.105041858 -3.733097998
25 -0.490714341 -0.105041858
26 -2.457174786 -0.490714341
27 0.002862262 -2.457174786
28 -1.540574800 0.002862262
29 -2.030615255 -1.540574800
30 7.200383000 -2.030615255
31 -1.107525711 7.200383000
32 -5.565736347 -1.107525711
33 0.827318374 -5.565736347
34 0.607835731 0.827318374
35 0.282515485 0.607835731
36 3.252829771 0.282515485
37 2.613320646 3.252829771
38 6.235355970 2.613320646
39 3.616627163 6.235355970
40 1.832031098 3.616627163
41 0.939980135 1.832031098
42 1.005932161 0.939980135
43 0.850704106 1.005932161
44 2.774387270 0.850704106
45 -3.650802594 2.774387270
46 -0.706810160 -3.650802594
47 -3.377079716 -0.706810160
48 -2.863722988 -3.377079716
49 -3.573125462 -2.863722988
50 -2.311668241 -3.573125462
51 -0.322469968 -2.311668241
52 -4.281113406 -0.322469968
53 1.557855816 -4.281113406
54 -2.763432878 1.557855816
55 1.397132327 -2.763432878
56 4.099615099 1.397132327
57 3.189270718 4.099615099
58 NA 3.189270718
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -1.891414074 -0.389998991
[2,] -0.643598441 -1.891414074
[3,] 1.427259050 -0.643598441
[4,] 0.936508116 1.427259050
[5,] 2.828172527 0.936508116
[6,] 6.581913218 2.828172527
[7,] 1.882652063 6.581913218
[8,] 1.786283920 1.882652063
[9,] -0.428332015 1.786283920
[10,] 2.025278451 -0.428332015
[11,] -0.029471374 2.025278451
[12,] -2.376073504 -0.029471374
[13,] 0.215528546 -2.376073504
[14,] -5.376151380 0.215528546
[15,] -1.709777753 -5.376151380
[16,] -1.501448029 -1.709777753
[17,] -0.359094720 -1.501448029
[18,] -7.069319050 -0.359094720
[19,] 3.918355672 -7.069319050
[20,] -2.268128703 3.918355672
[21,] 0.029875362 -2.268128703
[22,] 1.005730484 0.029875362
[23,] -3.733097998 1.005730484
[24,] -0.105041858 -3.733097998
[25,] -0.490714341 -0.105041858
[26,] -2.457174786 -0.490714341
[27,] 0.002862262 -2.457174786
[28,] -1.540574800 0.002862262
[29,] -2.030615255 -1.540574800
[30,] 7.200383000 -2.030615255
[31,] -1.107525711 7.200383000
[32,] -5.565736347 -1.107525711
[33,] 0.827318374 -5.565736347
[34,] 0.607835731 0.827318374
[35,] 0.282515485 0.607835731
[36,] 3.252829771 0.282515485
[37,] 2.613320646 3.252829771
[38,] 6.235355970 2.613320646
[39,] 3.616627163 6.235355970
[40,] 1.832031098 3.616627163
[41,] 0.939980135 1.832031098
[42,] 1.005932161 0.939980135
[43,] 0.850704106 1.005932161
[44,] 2.774387270 0.850704106
[45,] -3.650802594 2.774387270
[46,] -0.706810160 -3.650802594
[47,] -3.377079716 -0.706810160
[48,] -2.863722988 -3.377079716
[49,] -3.573125462 -2.863722988
[50,] -2.311668241 -3.573125462
[51,] -0.322469968 -2.311668241
[52,] -4.281113406 -0.322469968
[53,] 1.557855816 -4.281113406
[54,] -2.763432878 1.557855816
[55,] 1.397132327 -2.763432878
[56,] 4.099615099 1.397132327
[57,] 3.189270718 4.099615099
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -1.891414074 -0.389998991
2 -0.643598441 -1.891414074
3 1.427259050 -0.643598441
4 0.936508116 1.427259050
5 2.828172527 0.936508116
6 6.581913218 2.828172527
7 1.882652063 6.581913218
8 1.786283920 1.882652063
9 -0.428332015 1.786283920
10 2.025278451 -0.428332015
11 -0.029471374 2.025278451
12 -2.376073504 -0.029471374
13 0.215528546 -2.376073504
14 -5.376151380 0.215528546
15 -1.709777753 -5.376151380
16 -1.501448029 -1.709777753
17 -0.359094720 -1.501448029
18 -7.069319050 -0.359094720
19 3.918355672 -7.069319050
20 -2.268128703 3.918355672
21 0.029875362 -2.268128703
22 1.005730484 0.029875362
23 -3.733097998 1.005730484
24 -0.105041858 -3.733097998
25 -0.490714341 -0.105041858
26 -2.457174786 -0.490714341
27 0.002862262 -2.457174786
28 -1.540574800 0.002862262
29 -2.030615255 -1.540574800
30 7.200383000 -2.030615255
31 -1.107525711 7.200383000
32 -5.565736347 -1.107525711
33 0.827318374 -5.565736347
34 0.607835731 0.827318374
35 0.282515485 0.607835731
36 3.252829771 0.282515485
37 2.613320646 3.252829771
38 6.235355970 2.613320646
39 3.616627163 6.235355970
40 1.832031098 3.616627163
41 0.939980135 1.832031098
42 1.005932161 0.939980135
43 0.850704106 1.005932161
44 2.774387270 0.850704106
45 -3.650802594 2.774387270
46 -0.706810160 -3.650802594
47 -3.377079716 -0.706810160
48 -2.863722988 -3.377079716
49 -3.573125462 -2.863722988
50 -2.311668241 -3.573125462
51 -0.322469968 -2.311668241
52 -4.281113406 -0.322469968
53 1.557855816 -4.281113406
54 -2.763432878 1.557855816
55 1.397132327 -2.763432878
56 4.099615099 1.397132327
57 3.189270718 4.099615099
> 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/wessaorg/rcomp/tmp/7lhmm1321988471.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/wessaorg/rcomp/tmp/8yqoe1321988471.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/wessaorg/rcomp/tmp/9dwtg1321988471.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')
Warning messages:
1: In sqrt(crit * p * (1 - hh)/hh) : NaNs produced
2: In sqrt(crit * p * (1 - hh)/hh) : NaNs produced
> par(opar)
> dev.off()
null device
1
> if (n > n25) {
+ postscript(file="/var/wessaorg/rcomp/tmp/103s061321988471.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/wessaorg/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/wessaorg/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/wessaorg/rcomp/tmp/11j2fb1321988471.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/wessaorg/rcomp/tmp/12tywl1321988471.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/wessaorg/rcomp/tmp/13665l1321988471.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/wessaorg/rcomp/tmp/14kyo51321988471.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/wessaorg/rcomp/tmp/15gbre1321988471.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/wessaorg/rcomp/tmp/162gtf1321988471.tab")
+ }
>
> try(system("convert tmp/1gm8e1321988471.ps tmp/1gm8e1321988471.png",intern=TRUE))
character(0)
> try(system("convert tmp/227yk1321988471.ps tmp/227yk1321988471.png",intern=TRUE))
character(0)
> try(system("convert tmp/32q291321988471.ps tmp/32q291321988471.png",intern=TRUE))
character(0)
> try(system("convert tmp/4sqjk1321988471.ps tmp/4sqjk1321988471.png",intern=TRUE))
character(0)
> try(system("convert tmp/5s7ak1321988471.ps tmp/5s7ak1321988471.png",intern=TRUE))
character(0)
> try(system("convert tmp/6fet31321988471.ps tmp/6fet31321988471.png",intern=TRUE))
character(0)
> try(system("convert tmp/7lhmm1321988471.ps tmp/7lhmm1321988471.png",intern=TRUE))
character(0)
> try(system("convert tmp/8yqoe1321988471.ps tmp/8yqoe1321988471.png",intern=TRUE))
character(0)
> try(system("convert tmp/9dwtg1321988471.ps tmp/9dwtg1321988471.png",intern=TRUE))
character(0)
> try(system("convert tmp/103s061321988471.ps tmp/103s061321988471.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.327 0.530 3.901