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(286602,283042,276687,277915,277128,277103,275037,270150,267140,264993,287259,291186,292300,288186,281477,282656,280190,280408,276836,275216,274352,271311,289802,290726,292300,278506,269826,265861,269034,264176,255198,253353,246057,235372,258556,260993,254663,250643,243422,247105,248541,245039,237080,237085,225554,226839,247934,248333,246969,245098,246263,255765,264319,268347,273046,273963,267430,271993,292710,295881,294563),dim=c(1,61),dimnames=list(c('HPC'),1:61))
> y <- array(NA,dim=c(1,61),dimnames=list(c('HPC'),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 = '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
HPC M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 286602 1 0 0 0 0 0 0 0 0 0 0 1
2 283042 0 1 0 0 0 0 0 0 0 0 0 2
3 276687 0 0 1 0 0 0 0 0 0 0 0 3
4 277915 0 0 0 1 0 0 0 0 0 0 0 4
5 277128 0 0 0 0 1 0 0 0 0 0 0 5
6 277103 0 0 0 0 0 1 0 0 0 0 0 6
7 275037 0 0 0 0 0 0 1 0 0 0 0 7
8 270150 0 0 0 0 0 0 0 1 0 0 0 8
9 267140 0 0 0 0 0 0 0 0 1 0 0 9
10 264993 0 0 0 0 0 0 0 0 0 1 0 10
11 287259 0 0 0 0 0 0 0 0 0 0 1 11
12 291186 0 0 0 0 0 0 0 0 0 0 0 12
13 292300 1 0 0 0 0 0 0 0 0 0 0 13
14 288186 0 1 0 0 0 0 0 0 0 0 0 14
15 281477 0 0 1 0 0 0 0 0 0 0 0 15
16 282656 0 0 0 1 0 0 0 0 0 0 0 16
17 280190 0 0 0 0 1 0 0 0 0 0 0 17
18 280408 0 0 0 0 0 1 0 0 0 0 0 18
19 276836 0 0 0 0 0 0 1 0 0 0 0 19
20 275216 0 0 0 0 0 0 0 1 0 0 0 20
21 274352 0 0 0 0 0 0 0 0 1 0 0 21
22 271311 0 0 0 0 0 0 0 0 0 1 0 22
23 289802 0 0 0 0 0 0 0 0 0 0 1 23
24 290726 0 0 0 0 0 0 0 0 0 0 0 24
25 292300 1 0 0 0 0 0 0 0 0 0 0 25
26 278506 0 1 0 0 0 0 0 0 0 0 0 26
27 269826 0 0 1 0 0 0 0 0 0 0 0 27
28 265861 0 0 0 1 0 0 0 0 0 0 0 28
29 269034 0 0 0 0 1 0 0 0 0 0 0 29
30 264176 0 0 0 0 0 1 0 0 0 0 0 30
31 255198 0 0 0 0 0 0 1 0 0 0 0 31
32 253353 0 0 0 0 0 0 0 1 0 0 0 32
33 246057 0 0 0 0 0 0 0 0 1 0 0 33
34 235372 0 0 0 0 0 0 0 0 0 1 0 34
35 258556 0 0 0 0 0 0 0 0 0 0 1 35
36 260993 0 0 0 0 0 0 0 0 0 0 0 36
37 254663 1 0 0 0 0 0 0 0 0 0 0 37
38 250643 0 1 0 0 0 0 0 0 0 0 0 38
39 243422 0 0 1 0 0 0 0 0 0 0 0 39
40 247105 0 0 0 1 0 0 0 0 0 0 0 40
41 248541 0 0 0 0 1 0 0 0 0 0 0 41
42 245039 0 0 0 0 0 1 0 0 0 0 0 42
43 237080 0 0 0 0 0 0 1 0 0 0 0 43
44 237085 0 0 0 0 0 0 0 1 0 0 0 44
45 225554 0 0 0 0 0 0 0 0 1 0 0 45
46 226839 0 0 0 0 0 0 0 0 0 1 0 46
47 247934 0 0 0 0 0 0 0 0 0 0 1 47
48 248333 0 0 0 0 0 0 0 0 0 0 0 48
49 246969 1 0 0 0 0 0 0 0 0 0 0 49
50 245098 0 1 0 0 0 0 0 0 0 0 0 50
51 246263 0 0 1 0 0 0 0 0 0 0 0 51
52 255765 0 0 0 1 0 0 0 0 0 0 0 52
53 264319 0 0 0 0 1 0 0 0 0 0 0 53
54 268347 0 0 0 0 0 1 0 0 0 0 0 54
55 273046 0 0 0 0 0 0 1 0 0 0 0 55
56 273963 0 0 0 0 0 0 0 1 0 0 0 56
57 267430 0 0 0 0 0 0 0 0 1 0 0 57
58 271993 0 0 0 0 0 0 0 0 0 1 0 58
59 292710 0 0 0 0 0 0 0 0 0 0 1 59
60 295881 0 0 0 0 0 0 0 0 0 0 0 60
61 294563 1 0 0 0 0 0 0 0 0 0 0 61
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) M1 M2 M3 M4 M5
293568.4 -1766.6 -12813.4 -17925.0 -15151.1 -12720.6
M6 M7 M8 M9 M10 M11
-13100.0 -16226.7 -17264.2 -22662.6 -24219.1 -2620.1
t
-448.5
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-25171.1 -13233.9 834.5 9168.3 30117.4
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 293568.4 8513.6 34.482 < 2e-16 ***
M1 -1766.6 9928.8 -0.178 0.859528
M2 -12813.4 10421.3 -1.230 0.224863
M3 -17925.0 10408.0 -1.722 0.091467 .
M4 -15151.1 10396.1 -1.457 0.151522
M5 -12720.6 10385.6 -1.225 0.226613
M6 -13100.0 10376.4 -1.262 0.212878
M7 -16226.7 10368.7 -1.565 0.124159
M8 -17264.2 10362.4 -1.666 0.102217
M9 -22662.6 10357.4 -2.188 0.033566 *
M10 -24219.1 10353.9 -2.339 0.023539 *
M11 -2620.1 10351.8 -0.253 0.801270
t -448.5 120.8 -3.713 0.000533 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 16370 on 48 degrees of freedom
Multiple R-squared: 0.3523, Adjusted R-squared: 0.1904
F-statistic: 2.176 on 12 and 48 DF, p-value: 0.02862
> 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,] 6.016580e-06 1.203316e-05 0.999993983
[2,] 3.494089e-06 6.988178e-06 0.999996506
[3,] 2.545627e-07 5.091254e-07 0.999999745
[4,] 8.071250e-08 1.614250e-07 0.999999919
[5,] 4.783297e-09 9.566594e-09 0.999999995
[6,] 1.780231e-09 3.560462e-09 0.999999998
[7,] 2.068806e-10 4.137611e-10 1.000000000
[8,] 3.249367e-11 6.498734e-11 1.000000000
[9,] 6.678331e-11 1.335666e-10 1.000000000
[10,] 3.275842e-11 6.551684e-11 1.000000000
[11,] 3.513347e-08 7.026693e-08 0.999999965
[12,] 6.299526e-07 1.259905e-06 0.999999370
[13,] 9.995006e-06 1.999001e-05 0.999990005
[14,] 1.504443e-05 3.008886e-05 0.999984956
[15,] 4.587784e-05 9.175568e-05 0.999954122
[16,] 2.448968e-04 4.897937e-04 0.999755103
[17,] 5.352226e-04 1.070445e-03 0.999464777
[18,] 2.390384e-03 4.780768e-03 0.997609616
[19,] 9.891508e-03 1.978302e-02 0.990108492
[20,] 2.041223e-02 4.082445e-02 0.979587773
[21,] 4.293909e-02 8.587818e-02 0.957060910
[22,] 8.966125e-02 1.793225e-01 0.910338749
[23,] 2.410671e-01 4.821342e-01 0.758932885
[24,] 4.621065e-01 9.242131e-01 0.537893473
[25,] 7.223856e-01 5.552288e-01 0.277614383
[26,] 9.171755e-01 1.656489e-01 0.082824474
[27,] 9.922344e-01 1.553118e-02 0.007765591
[28,] 9.930024e-01 1.399530e-02 0.006997650
[29,] 9.979896e-01 4.020790e-03 0.002010395
[30,] 9.980006e-01 3.998794e-03 0.001999397
> postscript(file="/var/www/html/freestat/rcomp/tmp/17eqp1291924233.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/2i6qa1291924233.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/3i6qa1291924233.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/4i6qa1291924233.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/5i6qa1291924233.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 = 61
Frequency = 1
1 2 3 4 5 6
-4751.3725 3183.9020 2388.9020 1291.5020 -1477.4980 -674.6980
7 8 9 10 11 12
834.5020 -2566.4980 270.3020 128.3020 1243.7020 2999.1020
13 14 15 16 17 18
6328.1765 13709.4510 12560.4510 11414.0510 6966.0510 8011.8510
19 20 21 22 23 24
8015.0510 7881.0510 12863.8510 11827.8510 9168.2510 7920.6510
25 26 27 28 29 30
11709.7255 9411.0000 6291.0000 0.6000 1191.6000 -2838.6000
31 32 33 34 35 36
-8241.4000 -8600.4000 -10049.6000 -18729.6000 -16696.2000 -16430.8000
37 38 39 40 41 42
-20545.7255 -13070.4510 -14731.4510 -13373.8510 -13919.8510 -16594.0510
43 44 45 46 47 48
-20977.8510 -19486.8510 -25171.0510 -21881.0510 -21936.6510 -23709.2510
49 50 51 52 53 54
-22858.1765 -13233.9020 -6508.9020 667.6980 7239.6980 12095.4980
55 56 57 58 59 60
20369.6980 22772.6980 22086.4980 28654.4980 28220.8980 29220.2980
61
30117.3725
> postscript(file="/var/www/html/freestat/rcomp/tmp/6af7v1291924233.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 = 61
Frequency = 1
lag(myerror, k = 1) myerror
0 -4751.3725 NA
1 3183.9020 -4751.3725
2 2388.9020 3183.9020
3 1291.5020 2388.9020
4 -1477.4980 1291.5020
5 -674.6980 -1477.4980
6 834.5020 -674.6980
7 -2566.4980 834.5020
8 270.3020 -2566.4980
9 128.3020 270.3020
10 1243.7020 128.3020
11 2999.1020 1243.7020
12 6328.1765 2999.1020
13 13709.4510 6328.1765
14 12560.4510 13709.4510
15 11414.0510 12560.4510
16 6966.0510 11414.0510
17 8011.8510 6966.0510
18 8015.0510 8011.8510
19 7881.0510 8015.0510
20 12863.8510 7881.0510
21 11827.8510 12863.8510
22 9168.2510 11827.8510
23 7920.6510 9168.2510
24 11709.7255 7920.6510
25 9411.0000 11709.7255
26 6291.0000 9411.0000
27 0.6000 6291.0000
28 1191.6000 0.6000
29 -2838.6000 1191.6000
30 -8241.4000 -2838.6000
31 -8600.4000 -8241.4000
32 -10049.6000 -8600.4000
33 -18729.6000 -10049.6000
34 -16696.2000 -18729.6000
35 -16430.8000 -16696.2000
36 -20545.7255 -16430.8000
37 -13070.4510 -20545.7255
38 -14731.4510 -13070.4510
39 -13373.8510 -14731.4510
40 -13919.8510 -13373.8510
41 -16594.0510 -13919.8510
42 -20977.8510 -16594.0510
43 -19486.8510 -20977.8510
44 -25171.0510 -19486.8510
45 -21881.0510 -25171.0510
46 -21936.6510 -21881.0510
47 -23709.2510 -21936.6510
48 -22858.1765 -23709.2510
49 -13233.9020 -22858.1765
50 -6508.9020 -13233.9020
51 667.6980 -6508.9020
52 7239.6980 667.6980
53 12095.4980 7239.6980
54 20369.6980 12095.4980
55 22772.6980 20369.6980
56 22086.4980 22772.6980
57 28654.4980 22086.4980
58 28220.8980 28654.4980
59 29220.2980 28220.8980
60 30117.3725 29220.2980
61 NA 30117.3725
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 3183.9020 -4751.3725
[2,] 2388.9020 3183.9020
[3,] 1291.5020 2388.9020
[4,] -1477.4980 1291.5020
[5,] -674.6980 -1477.4980
[6,] 834.5020 -674.6980
[7,] -2566.4980 834.5020
[8,] 270.3020 -2566.4980
[9,] 128.3020 270.3020
[10,] 1243.7020 128.3020
[11,] 2999.1020 1243.7020
[12,] 6328.1765 2999.1020
[13,] 13709.4510 6328.1765
[14,] 12560.4510 13709.4510
[15,] 11414.0510 12560.4510
[16,] 6966.0510 11414.0510
[17,] 8011.8510 6966.0510
[18,] 8015.0510 8011.8510
[19,] 7881.0510 8015.0510
[20,] 12863.8510 7881.0510
[21,] 11827.8510 12863.8510
[22,] 9168.2510 11827.8510
[23,] 7920.6510 9168.2510
[24,] 11709.7255 7920.6510
[25,] 9411.0000 11709.7255
[26,] 6291.0000 9411.0000
[27,] 0.6000 6291.0000
[28,] 1191.6000 0.6000
[29,] -2838.6000 1191.6000
[30,] -8241.4000 -2838.6000
[31,] -8600.4000 -8241.4000
[32,] -10049.6000 -8600.4000
[33,] -18729.6000 -10049.6000
[34,] -16696.2000 -18729.6000
[35,] -16430.8000 -16696.2000
[36,] -20545.7255 -16430.8000
[37,] -13070.4510 -20545.7255
[38,] -14731.4510 -13070.4510
[39,] -13373.8510 -14731.4510
[40,] -13919.8510 -13373.8510
[41,] -16594.0510 -13919.8510
[42,] -20977.8510 -16594.0510
[43,] -19486.8510 -20977.8510
[44,] -25171.0510 -19486.8510
[45,] -21881.0510 -25171.0510
[46,] -21936.6510 -21881.0510
[47,] -23709.2510 -21936.6510
[48,] -22858.1765 -23709.2510
[49,] -13233.9020 -22858.1765
[50,] -6508.9020 -13233.9020
[51,] 667.6980 -6508.9020
[52,] 7239.6980 667.6980
[53,] 12095.4980 7239.6980
[54,] 20369.6980 12095.4980
[55,] 22772.6980 20369.6980
[56,] 22086.4980 22772.6980
[57,] 28654.4980 22086.4980
[58,] 28220.8980 28654.4980
[59,] 29220.2980 28220.8980
[60,] 30117.3725 29220.2980
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 3183.9020 -4751.3725
2 2388.9020 3183.9020
3 1291.5020 2388.9020
4 -1477.4980 1291.5020
5 -674.6980 -1477.4980
6 834.5020 -674.6980
7 -2566.4980 834.5020
8 270.3020 -2566.4980
9 128.3020 270.3020
10 1243.7020 128.3020
11 2999.1020 1243.7020
12 6328.1765 2999.1020
13 13709.4510 6328.1765
14 12560.4510 13709.4510
15 11414.0510 12560.4510
16 6966.0510 11414.0510
17 8011.8510 6966.0510
18 8015.0510 8011.8510
19 7881.0510 8015.0510
20 12863.8510 7881.0510
21 11827.8510 12863.8510
22 9168.2510 11827.8510
23 7920.6510 9168.2510
24 11709.7255 7920.6510
25 9411.0000 11709.7255
26 6291.0000 9411.0000
27 0.6000 6291.0000
28 1191.6000 0.6000
29 -2838.6000 1191.6000
30 -8241.4000 -2838.6000
31 -8600.4000 -8241.4000
32 -10049.6000 -8600.4000
33 -18729.6000 -10049.6000
34 -16696.2000 -18729.6000
35 -16430.8000 -16696.2000
36 -20545.7255 -16430.8000
37 -13070.4510 -20545.7255
38 -14731.4510 -13070.4510
39 -13373.8510 -14731.4510
40 -13919.8510 -13373.8510
41 -16594.0510 -13919.8510
42 -20977.8510 -16594.0510
43 -19486.8510 -20977.8510
44 -25171.0510 -19486.8510
45 -21881.0510 -25171.0510
46 -21936.6510 -21881.0510
47 -23709.2510 -21936.6510
48 -22858.1765 -23709.2510
49 -13233.9020 -22858.1765
50 -6508.9020 -13233.9020
51 667.6980 -6508.9020
52 7239.6980 667.6980
53 12095.4980 7239.6980
54 20369.6980 12095.4980
55 22772.6980 20369.6980
56 22086.4980 22772.6980
57 28654.4980 22086.4980
58 28220.8980 28654.4980
59 29220.2980 28220.8980
60 30117.3725 29220.2980
> 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/7l6of1291924233.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/8l6of1291924233.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/9l6of1291924233.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/10wx501291924233.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/11hgm61291924233.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/122y2c1291924233.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/13rzh61291924233.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/14pbnf1291924233.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/155rff1291924233.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/1611d51291924233.tab")
+ }
>
> try(system("convert tmp/17eqp1291924233.ps tmp/17eqp1291924233.png",intern=TRUE))
character(0)
> try(system("convert tmp/2i6qa1291924233.ps tmp/2i6qa1291924233.png",intern=TRUE))
character(0)
> try(system("convert tmp/3i6qa1291924233.ps tmp/3i6qa1291924233.png",intern=TRUE))
character(0)
> try(system("convert tmp/4i6qa1291924233.ps tmp/4i6qa1291924233.png",intern=TRUE))
character(0)
> try(system("convert tmp/5i6qa1291924233.ps tmp/5i6qa1291924233.png",intern=TRUE))
character(0)
> try(system("convert tmp/6af7v1291924233.ps tmp/6af7v1291924233.png",intern=TRUE))
character(0)
> try(system("convert tmp/7l6of1291924233.ps tmp/7l6of1291924233.png",intern=TRUE))
character(0)
> try(system("convert tmp/8l6of1291924233.ps tmp/8l6of1291924233.png",intern=TRUE))
character(0)
> try(system("convert tmp/9l6of1291924233.ps tmp/9l6of1291924233.png",intern=TRUE))
character(0)
> try(system("convert tmp/10wx501291924233.ps tmp/10wx501291924233.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.799 2.462 4.115