R version 2.8.0 (2008-10-20)
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.
Natural language support but running in an English locale
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(109.99
+ ,89
+ ,103.88
+ ,103.77
+ ,112.01
+ ,86.4
+ ,103.91
+ ,103.88
+ ,111.96
+ ,84.5
+ ,103.91
+ ,103.91
+ ,111.41
+ ,82.7
+ ,103.92
+ ,103.91
+ ,112.11
+ ,80.8
+ ,104.05
+ ,103.92
+ ,111.67
+ ,81.8
+ ,104.23
+ ,104.05
+ ,111.95
+ ,81.8
+ ,104.30
+ ,104.23
+ ,112.31
+ ,82.9
+ ,104.31
+ ,104.30
+ ,113.26
+ ,83.8
+ ,104.31
+ ,104.31
+ ,113.5
+ ,86.2
+ ,104.34
+ ,104.31
+ ,114.43
+ ,86.1
+ ,104.55
+ ,104.34
+ ,115.02
+ ,86.2
+ ,104.65
+ ,104.55
+ ,115.1
+ ,88.8
+ ,104.73
+ ,104.65
+ ,117.11
+ ,89.6
+ ,104.75
+ ,104.73
+ ,117.52
+ ,87.8
+ ,104.75
+ ,104.75
+ ,116.1
+ ,88.3
+ ,104.76
+ ,104.75
+ ,116.39
+ ,88.6
+ ,104.94
+ ,104.76
+ ,116.01
+ ,91
+ ,105.29
+ ,104.94
+ ,116.74
+ ,91.5
+ ,105.38
+ ,105.29
+ ,116.68
+ ,95.4
+ ,105.43
+ ,105.38
+ ,117.45
+ ,98.7
+ ,105.43
+ ,105.43
+ ,117.8
+ ,99.9
+ ,105.42
+ ,105.43
+ ,119.37
+ ,98.6
+ ,105.52
+ ,105.42
+ ,118.9
+ ,100.3
+ ,105.69
+ ,105.52
+ ,119.05
+ ,100.2
+ ,105.72
+ ,105.69
+ ,120.46
+ ,100.4
+ ,105.74
+ ,105.72
+ ,120.99
+ ,101.4
+ ,105.74
+ ,105.74
+ ,119.86
+ ,103
+ ,105.74
+ ,105.74
+ ,120.18
+ ,109.1
+ ,105.95
+ ,105.74
+ ,119.81
+ ,111.4
+ ,106.17
+ ,105.95
+ ,120.15
+ ,114.1
+ ,106.34
+ ,106.17
+ ,119.8
+ ,121.8
+ ,106.37
+ ,106.34
+ ,120.27
+ ,127.6
+ ,106.37
+ ,106.37
+ ,120.71
+ ,129.9
+ ,106.36
+ ,106.37
+ ,121.87
+ ,128
+ ,106.44
+ ,106.36
+ ,121.87
+ ,123.5
+ ,106.29
+ ,106.44
+ ,121.92
+ ,124
+ ,106.23
+ ,106.29
+ ,123.72
+ ,127.4
+ ,106.23
+ ,106.23
+ ,124.38
+ ,127.6
+ ,106.23
+ ,106.23
+ ,123.21
+ ,128.4
+ ,106.34
+ ,106.23
+ ,123.17
+ ,131.4
+ ,106.44
+ ,106.34
+ ,122.95
+ ,135.1
+ ,106.44
+ ,106.44
+ ,123.46
+ ,134
+ ,106.48
+ ,106.44
+ ,123.24
+ ,144.5
+ ,106.50
+ ,106.48
+ ,123.86
+ ,147.3
+ ,106.57
+ ,106.50
+ ,124.28
+ ,150.9
+ ,106.40
+ ,106.57
+ ,124.78
+ ,148.7
+ ,106.37
+ ,106.40
+ ,125.19
+ ,141.4
+ ,106.25
+ ,106.37
+ ,125.46
+ ,138.9
+ ,106.21
+ ,106.25
+ ,127.6
+ ,139.8
+ ,106.21
+ ,106.21
+ ,127.8
+ ,145.6
+ ,106.24
+ ,106.21
+ ,126.63
+ ,147.9
+ ,106.19
+ ,106.24
+ ,127.06
+ ,148.5
+ ,106.08
+ ,106.19
+ ,126.77
+ ,151.1
+ ,106.13
+ ,106.08)
+ ,dim=c(4
+ ,54)
+ ,dimnames=list(c('Y'
+ ,'X'
+ ,'Y1'
+ ,'Y2')
+ ,1:54))
> y <- array(NA,dim=c(4,54),dimnames=list(c('Y','X','Y1','Y2'),1:54))
> for (i in 1:dim(x)[1])
+ {
+ for (j in 1:dim(x)[2])
+ {
+ y[i,j] <- as.numeric(x[i,j])
+ }
+ }
> par3 = 'Linear Trend'
> par2 = 'Do not include Seasonal Dummies'
> par1 = '1'
> #'GNU S' R Code compiled by R2WASP v. 1.0.44 ()
> #Author: Prof. Dr. P. Wessa
> #To cite this work: AUTHOR(S), (YEAR), YOUR SOFTWARE TITLE (vNUMBER) in Free Statistics Software (v$_version), Office for Research Development and Education, URL http://www.wessa.net/rwasp_YOURPAGE.wasp/
> #Source of accompanying publication: Office for Research, Development, and Education
> #Technical description: Write here your technical program description (don't use hard returns!)
> library(lattice)
> library(lmtest)
Loading required package: zoo
Attaching package: 'zoo'
The following object(s) are masked from package:base :
as.Date.numeric
> n25 <- 25 #minimum number of obs. for Goldfeld-Quandt test
> par1 <- as.numeric(par1)
> x <- t(y)
> k <- length(x[1,])
> n <- length(x[,1])
> x1 <- cbind(x[,par1], x[,1:k!=par1])
> mycolnames <- c(colnames(x)[par1], colnames(x)[1:k!=par1])
> colnames(x1) <- mycolnames #colnames(x)[par1]
> x <- x1
> if (par3 == 'First Differences'){
+ x2 <- array(0, dim=c(n-1,k), dimnames=list(1:(n-1), paste('(1-B)',colnames(x),sep='')))
+ for (i in 1:n-1) {
+ for (j in 1:k) {
+ x2[i,j] <- x[i+1,j] - x[i,j]
+ }
+ }
+ x <- x2
+ }
> if (par2 == 'Include Monthly Dummies'){
+ x2 <- array(0, dim=c(n,11), dimnames=list(1:n, paste('M', seq(1:11), sep ='')))
+ for (i in 1:11){
+ x2[seq(i,n,12),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> if (par2 == 'Include Quarterly Dummies'){
+ x2 <- array(0, dim=c(n,3), dimnames=list(1:n, paste('Q', seq(1:3), sep ='')))
+ for (i in 1:3){
+ x2[seq(i,n,4),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> k <- length(x[1,])
> if (par3 == 'Linear Trend'){
+ x <- cbind(x, c(1:n))
+ colnames(x)[k+1] <- 't'
+ }
> x
Y X Y1 Y2 t
1 109.99 89.0 103.88 103.77 1
2 112.01 86.4 103.91 103.88 2
3 111.96 84.5 103.91 103.91 3
4 111.41 82.7 103.92 103.91 4
5 112.11 80.8 104.05 103.92 5
6 111.67 81.8 104.23 104.05 6
7 111.95 81.8 104.30 104.23 7
8 112.31 82.9 104.31 104.30 8
9 113.26 83.8 104.31 104.31 9
10 113.50 86.2 104.34 104.31 10
11 114.43 86.1 104.55 104.34 11
12 115.02 86.2 104.65 104.55 12
13 115.10 88.8 104.73 104.65 13
14 117.11 89.6 104.75 104.73 14
15 117.52 87.8 104.75 104.75 15
16 116.10 88.3 104.76 104.75 16
17 116.39 88.6 104.94 104.76 17
18 116.01 91.0 105.29 104.94 18
19 116.74 91.5 105.38 105.29 19
20 116.68 95.4 105.43 105.38 20
21 117.45 98.7 105.43 105.43 21
22 117.80 99.9 105.42 105.43 22
23 119.37 98.6 105.52 105.42 23
24 118.90 100.3 105.69 105.52 24
25 119.05 100.2 105.72 105.69 25
26 120.46 100.4 105.74 105.72 26
27 120.99 101.4 105.74 105.74 27
28 119.86 103.0 105.74 105.74 28
29 120.18 109.1 105.95 105.74 29
30 119.81 111.4 106.17 105.95 30
31 120.15 114.1 106.34 106.17 31
32 119.80 121.8 106.37 106.34 32
33 120.27 127.6 106.37 106.37 33
34 120.71 129.9 106.36 106.37 34
35 121.87 128.0 106.44 106.36 35
36 121.87 123.5 106.29 106.44 36
37 121.92 124.0 106.23 106.29 37
38 123.72 127.4 106.23 106.23 38
39 124.38 127.6 106.23 106.23 39
40 123.21 128.4 106.34 106.23 40
41 123.17 131.4 106.44 106.34 41
42 122.95 135.1 106.44 106.44 42
43 123.46 134.0 106.48 106.44 43
44 123.24 144.5 106.50 106.48 44
45 123.86 147.3 106.57 106.50 45
46 124.28 150.9 106.40 106.57 46
47 124.78 148.7 106.37 106.40 47
48 125.19 141.4 106.25 106.37 48
49 125.46 138.9 106.21 106.25 49
50 127.60 139.8 106.21 106.21 50
51 127.80 145.6 106.24 106.21 51
52 126.63 147.9 106.19 106.24 52
53 127.06 148.5 106.08 106.19 53
54 126.77 151.1 106.13 106.08 54
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X Y1 Y2 t
118.42776 -0.08013 -1.64514 1.62948 0.42097
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-1.1933 -0.5256 -0.1643 0.2780 1.6410
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 118.42776 31.98194 3.703 0.000541 ***
X -0.08013 0.01740 -4.606 2.94e-05 ***
Y1 -1.64514 1.18102 -1.393 0.169915
Y2 1.62948 1.22591 1.329 0.189938
t 0.42097 0.03212 13.106 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.7301 on 49 degrees of freedom
Multiple R-squared: 0.9794, Adjusted R-squared: 0.9777
F-statistic: 581.3 on 4 and 49 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.4253518 0.8507036 0.5746482
[2,] 0.5339575 0.9320851 0.4660425
[3,] 0.4323912 0.8647823 0.5676088
[4,] 0.5463845 0.9072309 0.4536155
[5,] 0.5528277 0.8943446 0.4471723
[6,] 0.4403387 0.8806773 0.5596613
[7,] 0.6571969 0.6856062 0.3428031
[8,] 0.7758913 0.4482174 0.2241087
[9,] 0.7934929 0.4130141 0.2065071
[10,] 0.7319341 0.5361318 0.2680659
[11,] 0.7062169 0.5875662 0.2937831
[12,] 0.6981724 0.6036552 0.3018276
[13,] 0.7821333 0.4357335 0.2178667
[14,] 0.7796483 0.4407035 0.2203517
[15,] 0.7938383 0.4123235 0.2061617
[16,] 0.7638159 0.4723682 0.2361841
[17,] 0.6945099 0.6109801 0.3054901
[18,] 0.6496364 0.7007271 0.3503636
[19,] 0.6043121 0.7913759 0.3956879
[20,] 0.5938794 0.8122411 0.4061206
[21,] 0.6688225 0.6623550 0.3311775
[22,] 0.6162490 0.7675021 0.3837510
[23,] 0.6314047 0.7371906 0.3685953
[24,] 0.6045103 0.7909794 0.3954897
[25,] 0.6580112 0.6839776 0.3419888
[26,] 0.6238834 0.7522333 0.3761166
[27,] 0.5871192 0.8257616 0.4128808
[28,] 0.5097089 0.9805822 0.4902911
[29,] 0.4604353 0.9208705 0.5395647
[30,] 0.5545342 0.8909316 0.4454658
[31,] 0.4879935 0.9759870 0.5120065
[32,] 0.6090578 0.7818844 0.3909422
[33,] 0.5679445 0.8641111 0.4320555
[34,] 0.5307764 0.9384472 0.4692236
[35,] 0.4686186 0.9372372 0.5313814
[36,] 0.3967199 0.7934397 0.6032801
[37,] 0.2983147 0.5966294 0.7016853
[38,] 0.3608812 0.7217623 0.6391188
[39,] 0.2245884 0.4491768 0.7754116
> postscript(file="/var/www/html/freestat/rcomp/tmp/1dbu11293029226.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/freestat/rcomp/tmp/2dbu11293029226.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/freestat/rcomp/tmp/363b41293029226.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/freestat/rcomp/tmp/463b41293029226.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/freestat/rcomp/tmp/563b41293029226.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 = 54
Frequency = 1
1 2 3 4 5 6
0.07856148 1.33934943 0.66723596 -0.43152824 -0.10718442 -0.80373186
7 8 9 10 11 12
-1.12285274 -1.19329087 -0.60843832 -0.54773524 0.24987174 0.24923391
13 14 15 16 17 18
0.08527244 1.64095052 1.45314532 -0.33130971 -0.15841340 -0.48457369
19 20 21 22 23 24
-0.55773661 -0.79058257 -0.25858669 -0.24985035 0.97580969 0.33778945
25 26 27 28 29 30
-0.16885548 0.82021606 0.97678723 -0.44597129 0.28735355 -0.29957228
31 32 33 34 35 36
-0.24299595 -0.62459240 -0.15967083 0.02721338 0.76189005 -0.39681735
37 38 39 40 41 42
-0.58200957 1.16724278 1.42229607 0.07639492 -0.15890489 -0.66632934
43 44 45 46 47 48
-0.59964538 -0.43148414 0.07448849 -0.03173801 0.09865055 -0.64583613
49 50 51 52 53 54
-0.86741335 0.98891331 1.28207341 -0.25573227 -0.29811603 -0.53924033
> postscript(file="/var/www/html/freestat/rcomp/tmp/6zut71293029226.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 = 54
Frequency = 1
lag(myerror, k = 1) myerror
0 0.07856148 NA
1 1.33934943 0.07856148
2 0.66723596 1.33934943
3 -0.43152824 0.66723596
4 -0.10718442 -0.43152824
5 -0.80373186 -0.10718442
6 -1.12285274 -0.80373186
7 -1.19329087 -1.12285274
8 -0.60843832 -1.19329087
9 -0.54773524 -0.60843832
10 0.24987174 -0.54773524
11 0.24923391 0.24987174
12 0.08527244 0.24923391
13 1.64095052 0.08527244
14 1.45314532 1.64095052
15 -0.33130971 1.45314532
16 -0.15841340 -0.33130971
17 -0.48457369 -0.15841340
18 -0.55773661 -0.48457369
19 -0.79058257 -0.55773661
20 -0.25858669 -0.79058257
21 -0.24985035 -0.25858669
22 0.97580969 -0.24985035
23 0.33778945 0.97580969
24 -0.16885548 0.33778945
25 0.82021606 -0.16885548
26 0.97678723 0.82021606
27 -0.44597129 0.97678723
28 0.28735355 -0.44597129
29 -0.29957228 0.28735355
30 -0.24299595 -0.29957228
31 -0.62459240 -0.24299595
32 -0.15967083 -0.62459240
33 0.02721338 -0.15967083
34 0.76189005 0.02721338
35 -0.39681735 0.76189005
36 -0.58200957 -0.39681735
37 1.16724278 -0.58200957
38 1.42229607 1.16724278
39 0.07639492 1.42229607
40 -0.15890489 0.07639492
41 -0.66632934 -0.15890489
42 -0.59964538 -0.66632934
43 -0.43148414 -0.59964538
44 0.07448849 -0.43148414
45 -0.03173801 0.07448849
46 0.09865055 -0.03173801
47 -0.64583613 0.09865055
48 -0.86741335 -0.64583613
49 0.98891331 -0.86741335
50 1.28207341 0.98891331
51 -0.25573227 1.28207341
52 -0.29811603 -0.25573227
53 -0.53924033 -0.29811603
54 NA -0.53924033
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 1.33934943 0.07856148
[2,] 0.66723596 1.33934943
[3,] -0.43152824 0.66723596
[4,] -0.10718442 -0.43152824
[5,] -0.80373186 -0.10718442
[6,] -1.12285274 -0.80373186
[7,] -1.19329087 -1.12285274
[8,] -0.60843832 -1.19329087
[9,] -0.54773524 -0.60843832
[10,] 0.24987174 -0.54773524
[11,] 0.24923391 0.24987174
[12,] 0.08527244 0.24923391
[13,] 1.64095052 0.08527244
[14,] 1.45314532 1.64095052
[15,] -0.33130971 1.45314532
[16,] -0.15841340 -0.33130971
[17,] -0.48457369 -0.15841340
[18,] -0.55773661 -0.48457369
[19,] -0.79058257 -0.55773661
[20,] -0.25858669 -0.79058257
[21,] -0.24985035 -0.25858669
[22,] 0.97580969 -0.24985035
[23,] 0.33778945 0.97580969
[24,] -0.16885548 0.33778945
[25,] 0.82021606 -0.16885548
[26,] 0.97678723 0.82021606
[27,] -0.44597129 0.97678723
[28,] 0.28735355 -0.44597129
[29,] -0.29957228 0.28735355
[30,] -0.24299595 -0.29957228
[31,] -0.62459240 -0.24299595
[32,] -0.15967083 -0.62459240
[33,] 0.02721338 -0.15967083
[34,] 0.76189005 0.02721338
[35,] -0.39681735 0.76189005
[36,] -0.58200957 -0.39681735
[37,] 1.16724278 -0.58200957
[38,] 1.42229607 1.16724278
[39,] 0.07639492 1.42229607
[40,] -0.15890489 0.07639492
[41,] -0.66632934 -0.15890489
[42,] -0.59964538 -0.66632934
[43,] -0.43148414 -0.59964538
[44,] 0.07448849 -0.43148414
[45,] -0.03173801 0.07448849
[46,] 0.09865055 -0.03173801
[47,] -0.64583613 0.09865055
[48,] -0.86741335 -0.64583613
[49,] 0.98891331 -0.86741335
[50,] 1.28207341 0.98891331
[51,] -0.25573227 1.28207341
[52,] -0.29811603 -0.25573227
[53,] -0.53924033 -0.29811603
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 1.33934943 0.07856148
2 0.66723596 1.33934943
3 -0.43152824 0.66723596
4 -0.10718442 -0.43152824
5 -0.80373186 -0.10718442
6 -1.12285274 -0.80373186
7 -1.19329087 -1.12285274
8 -0.60843832 -1.19329087
9 -0.54773524 -0.60843832
10 0.24987174 -0.54773524
11 0.24923391 0.24987174
12 0.08527244 0.24923391
13 1.64095052 0.08527244
14 1.45314532 1.64095052
15 -0.33130971 1.45314532
16 -0.15841340 -0.33130971
17 -0.48457369 -0.15841340
18 -0.55773661 -0.48457369
19 -0.79058257 -0.55773661
20 -0.25858669 -0.79058257
21 -0.24985035 -0.25858669
22 0.97580969 -0.24985035
23 0.33778945 0.97580969
24 -0.16885548 0.33778945
25 0.82021606 -0.16885548
26 0.97678723 0.82021606
27 -0.44597129 0.97678723
28 0.28735355 -0.44597129
29 -0.29957228 0.28735355
30 -0.24299595 -0.29957228
31 -0.62459240 -0.24299595
32 -0.15967083 -0.62459240
33 0.02721338 -0.15967083
34 0.76189005 0.02721338
35 -0.39681735 0.76189005
36 -0.58200957 -0.39681735
37 1.16724278 -0.58200957
38 1.42229607 1.16724278
39 0.07639492 1.42229607
40 -0.15890489 0.07639492
41 -0.66632934 -0.15890489
42 -0.59964538 -0.66632934
43 -0.43148414 -0.59964538
44 0.07448849 -0.43148414
45 -0.03173801 0.07448849
46 0.09865055 -0.03173801
47 -0.64583613 0.09865055
48 -0.86741335 -0.64583613
49 0.98891331 -0.86741335
50 1.28207341 0.98891331
51 -0.25573227 1.28207341
52 -0.29811603 -0.25573227
53 -0.53924033 -0.29811603
> 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/freestat/rcomp/tmp/79lss1293029226.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/freestat/rcomp/tmp/89lss1293029226.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/freestat/rcomp/tmp/99lss1293029226.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')
> par(opar)
> dev.off()
null device
1
> if (n > n25) {
+ postscript(file="/var/www/html/freestat/rcomp/tmp/10kc9v1293029226.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/freestat/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/html/freestat/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/freestat/rcomp/tmp/11nd8j1293029226.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/freestat/rcomp/tmp/12159j1293029227.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/freestat/rcomp/tmp/13xf6s1293029227.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/freestat/rcomp/tmp/14iyng1293029227.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/freestat/rcomp/tmp/15my3m1293029227.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/freestat/rcomp/tmp/16phka1293029227.tab")
+ }
>
> try(system("convert tmp/1dbu11293029226.ps tmp/1dbu11293029226.png",intern=TRUE))
character(0)
> try(system("convert tmp/2dbu11293029226.ps tmp/2dbu11293029226.png",intern=TRUE))
character(0)
> try(system("convert tmp/363b41293029226.ps tmp/363b41293029226.png",intern=TRUE))
character(0)
> try(system("convert tmp/463b41293029226.ps tmp/463b41293029226.png",intern=TRUE))
character(0)
> try(system("convert tmp/563b41293029226.ps tmp/563b41293029226.png",intern=TRUE))
character(0)
> try(system("convert tmp/6zut71293029226.ps tmp/6zut71293029226.png",intern=TRUE))
character(0)
> try(system("convert tmp/79lss1293029226.ps tmp/79lss1293029226.png",intern=TRUE))
character(0)
> try(system("convert tmp/89lss1293029226.ps tmp/89lss1293029226.png",intern=TRUE))
character(0)
> try(system("convert tmp/99lss1293029226.ps tmp/99lss1293029226.png",intern=TRUE))
character(0)
> try(system("convert tmp/10kc9v1293029226.ps tmp/10kc9v1293029226.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.761 2.456 4.085