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(2253,14.9,2218,18.6,1855,19.1,2187,18.8,1852,18.2,1570,18,1851,19,1954,20.7,1828,21.2,2251,20.7,2277,19.6,2085,18.6,2282,18.7,2266,23.8,1878,24.9,2267,24.8,2069,23.8,1746,22.3,2299,21.7,2360,20.7,2214,19.7,2825,18.4,2355,17.4,2333,17,3016,18,2155,23.8,2172,25.5,2150,25.6,2533,23.7,2058,22,2160,21.3,2260,20.7,2498,20.4,2695,20.3,2799,20.4,2946,19.8,2930,19.5,2318,23.1,2540,23.5,2570,23.5,2669,22.9,2450,21.9,2842,21.5,3440,20.5,2678,20.2,2981,19.4,2260,19.2,2844,18.8,2546,18.8,2456,22.6,2295,23.3,2379,23,2479,21.4,2057,19.9,2280,18.8,2351,18.6,2276,18.4,2548,18.6,2311,19.9,2201,19.2),dim=c(2,60),dimnames=list(c('wngb','<25'),1:60))
> y <- array(NA,dim=c(2,60),dimnames=list(c('wngb','<25'),1:60))
> for (i in 1:dim(x)[1])
+ {
+ for (j in 1:dim(x)[2])
+ {
+ y[i,j] <- as.numeric(x[i,j])
+ }
+ }
> par3 = 'No Linear Trend'
> par2 = 'Include Monthly Dummies'
> par1 = '1'
> #'GNU S' R Code compiled by R2WASP v. 1.0.44 ()
> #Author: Prof. Dr. P. Wessa
> #To cite this work: AUTHOR(S), (YEAR), YOUR SOFTWARE TITLE (vNUMBER) in Free Statistics Software (v$_version), Office for Research Development and Education, URL http://www.wessa.net/rwasp_YOURPAGE.wasp/
> #Source of accompanying publication: Office for Research, Development, and Education
> #Technical description: Write here your technical program description (don't use hard returns!)
> library(lattice)
> library(lmtest)
Loading required package: zoo
Attaching package: 'zoo'
The following object(s) are masked from package:base :
as.Date.numeric
> n25 <- 25 #minimum number of obs. for Goldfeld-Quandt test
> par1 <- as.numeric(par1)
> x <- t(y)
> k <- length(x[1,])
> n <- length(x[,1])
> x1 <- cbind(x[,par1], x[,1:k!=par1])
> mycolnames <- c(colnames(x)[par1], colnames(x)[1:k!=par1])
> colnames(x1) <- mycolnames #colnames(x)[par1]
> x <- x1
> if (par3 == 'First Differences'){
+ x2 <- array(0, dim=c(n-1,k), dimnames=list(1:(n-1), paste('(1-B)',colnames(x),sep='')))
+ for (i in 1:n-1) {
+ for (j in 1:k) {
+ x2[i,j] <- x[i+1,j] - x[i,j]
+ }
+ }
+ x <- x2
+ }
> if (par2 == 'Include Monthly Dummies'){
+ x2 <- array(0, dim=c(n,11), dimnames=list(1:n, paste('M', seq(1:11), sep ='')))
+ for (i in 1:11){
+ x2[seq(i,n,12),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> if (par2 == 'Include Quarterly Dummies'){
+ x2 <- array(0, dim=c(n,3), dimnames=list(1:n, paste('Q', seq(1:3), sep ='')))
+ for (i in 1:3){
+ x2[seq(i,n,4),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> k <- length(x[1,])
> if (par3 == 'Linear Trend'){
+ x <- cbind(x, c(1:n))
+ colnames(x)[k+1] <- 't'
+ }
> x
wngb <25 M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 2253 14.9 1 0 0 0 0 0 0 0 0 0 0
2 2218 18.6 0 1 0 0 0 0 0 0 0 0 0
3 1855 19.1 0 0 1 0 0 0 0 0 0 0 0
4 2187 18.8 0 0 0 1 0 0 0 0 0 0 0
5 1852 18.2 0 0 0 0 1 0 0 0 0 0 0
6 1570 18.0 0 0 0 0 0 1 0 0 0 0 0
7 1851 19.0 0 0 0 0 0 0 1 0 0 0 0
8 1954 20.7 0 0 0 0 0 0 0 1 0 0 0
9 1828 21.2 0 0 0 0 0 0 0 0 1 0 0
10 2251 20.7 0 0 0 0 0 0 0 0 0 1 0
11 2277 19.6 0 0 0 0 0 0 0 0 0 0 1
12 2085 18.6 0 0 0 0 0 0 0 0 0 0 0
13 2282 18.7 1 0 0 0 0 0 0 0 0 0 0
14 2266 23.8 0 1 0 0 0 0 0 0 0 0 0
15 1878 24.9 0 0 1 0 0 0 0 0 0 0 0
16 2267 24.8 0 0 0 1 0 0 0 0 0 0 0
17 2069 23.8 0 0 0 0 1 0 0 0 0 0 0
18 1746 22.3 0 0 0 0 0 1 0 0 0 0 0
19 2299 21.7 0 0 0 0 0 0 1 0 0 0 0
20 2360 20.7 0 0 0 0 0 0 0 1 0 0 0
21 2214 19.7 0 0 0 0 0 0 0 0 1 0 0
22 2825 18.4 0 0 0 0 0 0 0 0 0 1 0
23 2355 17.4 0 0 0 0 0 0 0 0 0 0 1
24 2333 17.0 0 0 0 0 0 0 0 0 0 0 0
25 3016 18.0 1 0 0 0 0 0 0 0 0 0 0
26 2155 23.8 0 1 0 0 0 0 0 0 0 0 0
27 2172 25.5 0 0 1 0 0 0 0 0 0 0 0
28 2150 25.6 0 0 0 1 0 0 0 0 0 0 0
29 2533 23.7 0 0 0 0 1 0 0 0 0 0 0
30 2058 22.0 0 0 0 0 0 1 0 0 0 0 0
31 2160 21.3 0 0 0 0 0 0 1 0 0 0 0
32 2260 20.7 0 0 0 0 0 0 0 1 0 0 0
33 2498 20.4 0 0 0 0 0 0 0 0 1 0 0
34 2695 20.3 0 0 0 0 0 0 0 0 0 1 0
35 2799 20.4 0 0 0 0 0 0 0 0 0 0 1
36 2946 19.8 0 0 0 0 0 0 0 0 0 0 0
37 2930 19.5 1 0 0 0 0 0 0 0 0 0 0
38 2318 23.1 0 1 0 0 0 0 0 0 0 0 0
39 2540 23.5 0 0 1 0 0 0 0 0 0 0 0
40 2570 23.5 0 0 0 1 0 0 0 0 0 0 0
41 2669 22.9 0 0 0 0 1 0 0 0 0 0 0
42 2450 21.9 0 0 0 0 0 1 0 0 0 0 0
43 2842 21.5 0 0 0 0 0 0 1 0 0 0 0
44 3440 20.5 0 0 0 0 0 0 0 1 0 0 0
45 2678 20.2 0 0 0 0 0 0 0 0 1 0 0
46 2981 19.4 0 0 0 0 0 0 0 0 0 1 0
47 2260 19.2 0 0 0 0 0 0 0 0 0 0 1
48 2844 18.8 0 0 0 0 0 0 0 0 0 0 0
49 2546 18.8 1 0 0 0 0 0 0 0 0 0 0
50 2456 22.6 0 1 0 0 0 0 0 0 0 0 0
51 2295 23.3 0 0 1 0 0 0 0 0 0 0 0
52 2379 23.0 0 0 0 1 0 0 0 0 0 0 0
53 2479 21.4 0 0 0 0 1 0 0 0 0 0 0
54 2057 19.9 0 0 0 0 0 1 0 0 0 0 0
55 2280 18.8 0 0 0 0 0 0 1 0 0 0 0
56 2351 18.6 0 0 0 0 0 0 0 1 0 0 0
57 2276 18.4 0 0 0 0 0 0 0 0 1 0 0
58 2548 18.6 0 0 0 0 0 0 0 0 0 1 0
59 2311 19.9 0 0 0 0 0 0 0 0 0 0 1
60 2201 19.2 0 0 0 0 0 0 0 0 0 0 0
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) `<25` M1 M2 M3 M4
1588.51 47.82 157.07 -376.14 -552.82 -384.48
M5 M6 M7 M8 M9 M10
-320.16 -607.94 -280.52 -83.40 -245.17 139.94
M11
-111.05
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-541.00 -197.91 -45.14 181.16 954.57
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1588.51 514.80 3.086 0.00340 **
`<25` 47.82 26.45 1.808 0.07699 .
M1 157.07 205.48 0.764 0.44843
M2 -376.14 226.84 -1.658 0.10394
M3 -552.82 237.81 -2.325 0.02446 *
M4 -384.48 236.21 -1.628 0.11027
M5 -320.16 222.69 -1.438 0.15713
M6 -607.94 212.33 -2.863 0.00625 **
M7 -280.52 209.99 -1.336 0.18802
M8 -83.40 208.76 -0.399 0.69133
M9 -245.17 207.51 -1.181 0.24336
M10 139.94 205.73 0.680 0.49971
M11 -111.05 205.30 -0.541 0.59112
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 323.6 on 47 degrees of freedom
Multiple R-squared: 0.318, Adjusted R-squared: 0.1438
F-statistic: 1.826 on 12 and 47 DF, p-value: 0.07105
> 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.0002147164 0.0004294328 0.9997853
[2,] 0.0023271702 0.0046543403 0.9976728
[3,] 0.0010411482 0.0020822963 0.9989589
[4,] 0.0200808938 0.0401617875 0.9799191
[5,] 0.0492190508 0.0984381017 0.9507809
[6,] 0.0785731498 0.1571462995 0.9214269
[7,] 0.1837921306 0.3675842611 0.8162079
[8,] 0.1200079248 0.2400158497 0.8799921
[9,] 0.0947686163 0.1895372326 0.9052314
[10,] 0.2932179020 0.5864358040 0.7067821
[11,] 0.2312328742 0.4624657485 0.7687671
[12,] 0.2111582971 0.4223165941 0.7888417
[13,] 0.2137087961 0.4274175922 0.7862912
[14,] 0.2613361845 0.5226723691 0.7386638
[15,] 0.2631480962 0.5262961923 0.7368519
[16,] 0.2927766255 0.5855532510 0.7072234
[17,] 0.5745366410 0.8509267180 0.4254634
[18,] 0.6089646845 0.7820706310 0.3910353
[19,] 0.6335681847 0.7328636306 0.3664318
[20,] 0.6538034885 0.6923930230 0.3461965
[21,] 0.7110787570 0.5778424860 0.2889212
[22,] 0.6625284999 0.6749430002 0.3374715
[23,] 0.5972578650 0.8054842700 0.4027421
[24,] 0.5837997397 0.8324005207 0.4162003
[25,] 0.4921649642 0.9843299284 0.5078350
[26,] 0.4346448896 0.8692897793 0.5653551
[27,] 0.3801741617 0.7603483234 0.6198258
[28,] 0.3552330404 0.7104660807 0.6447670
[29,] 0.5289804490 0.9420391020 0.4710196
> postscript(file="/var/www/html/rcomp/tmp/1hrl11258735161.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/21xzq1258735161.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/3a3mq1258735161.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/47u061258735161.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/58qcr1258735161.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 = 60
Frequency = 1
1 2 3 4 5 6
-205.1124723 116.1619658 -94.0661963 83.9415163 -286.6816216 -271.3458350
7 8 9 10 11 12
-365.5818862 -540.9974879 -529.1411636 -467.3411636 -137.7461878 -392.9743499
13 14 15 16 17 18
-357.8308506 -84.5052888 -348.4258265 -122.9822390 -337.4771266 -300.9745263
19 20 21 22 23 24
-46.6975761 -134.9974879 -71.4102248 216.6462760 45.4591892 -68.4613485
25 26 27 28 29 30
409.6435875 -195.5052888 -83.1182020 -278.2387397 131.3049360 25.3716615
31 32 33 34 35 36
-166.5693257 -234.9974879 179.1153371 -4.2129132 345.9973115 410.6408990
37 38 39 40 41 42
251.9126487 0.9691494 380.5230498 242.1845747 305.5614367 422.1537240
43 44 45 46 47 48
505.8665491 954.5666373 368.6794623 324.8256501 -135.6179374 356.4615249
49 50 51 52 53 54
-98.6129132 162.8794623 145.0871750 75.0948876 187.2923755 124.7949758
55 56 57 58 59 60
72.9822390 -43.5741735 52.7565889 -69.9178492 -118.0923755 -305.6667255
> postscript(file="/var/www/html/rcomp/tmp/6gki61258735161.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 = 60
Frequency = 1
lag(myerror, k = 1) myerror
0 -205.1124723 NA
1 116.1619658 -205.1124723
2 -94.0661963 116.1619658
3 83.9415163 -94.0661963
4 -286.6816216 83.9415163
5 -271.3458350 -286.6816216
6 -365.5818862 -271.3458350
7 -540.9974879 -365.5818862
8 -529.1411636 -540.9974879
9 -467.3411636 -529.1411636
10 -137.7461878 -467.3411636
11 -392.9743499 -137.7461878
12 -357.8308506 -392.9743499
13 -84.5052888 -357.8308506
14 -348.4258265 -84.5052888
15 -122.9822390 -348.4258265
16 -337.4771266 -122.9822390
17 -300.9745263 -337.4771266
18 -46.6975761 -300.9745263
19 -134.9974879 -46.6975761
20 -71.4102248 -134.9974879
21 216.6462760 -71.4102248
22 45.4591892 216.6462760
23 -68.4613485 45.4591892
24 409.6435875 -68.4613485
25 -195.5052888 409.6435875
26 -83.1182020 -195.5052888
27 -278.2387397 -83.1182020
28 131.3049360 -278.2387397
29 25.3716615 131.3049360
30 -166.5693257 25.3716615
31 -234.9974879 -166.5693257
32 179.1153371 -234.9974879
33 -4.2129132 179.1153371
34 345.9973115 -4.2129132
35 410.6408990 345.9973115
36 251.9126487 410.6408990
37 0.9691494 251.9126487
38 380.5230498 0.9691494
39 242.1845747 380.5230498
40 305.5614367 242.1845747
41 422.1537240 305.5614367
42 505.8665491 422.1537240
43 954.5666373 505.8665491
44 368.6794623 954.5666373
45 324.8256501 368.6794623
46 -135.6179374 324.8256501
47 356.4615249 -135.6179374
48 -98.6129132 356.4615249
49 162.8794623 -98.6129132
50 145.0871750 162.8794623
51 75.0948876 145.0871750
52 187.2923755 75.0948876
53 124.7949758 187.2923755
54 72.9822390 124.7949758
55 -43.5741735 72.9822390
56 52.7565889 -43.5741735
57 -69.9178492 52.7565889
58 -118.0923755 -69.9178492
59 -305.6667255 -118.0923755
60 NA -305.6667255
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 116.1619658 -205.1124723
[2,] -94.0661963 116.1619658
[3,] 83.9415163 -94.0661963
[4,] -286.6816216 83.9415163
[5,] -271.3458350 -286.6816216
[6,] -365.5818862 -271.3458350
[7,] -540.9974879 -365.5818862
[8,] -529.1411636 -540.9974879
[9,] -467.3411636 -529.1411636
[10,] -137.7461878 -467.3411636
[11,] -392.9743499 -137.7461878
[12,] -357.8308506 -392.9743499
[13,] -84.5052888 -357.8308506
[14,] -348.4258265 -84.5052888
[15,] -122.9822390 -348.4258265
[16,] -337.4771266 -122.9822390
[17,] -300.9745263 -337.4771266
[18,] -46.6975761 -300.9745263
[19,] -134.9974879 -46.6975761
[20,] -71.4102248 -134.9974879
[21,] 216.6462760 -71.4102248
[22,] 45.4591892 216.6462760
[23,] -68.4613485 45.4591892
[24,] 409.6435875 -68.4613485
[25,] -195.5052888 409.6435875
[26,] -83.1182020 -195.5052888
[27,] -278.2387397 -83.1182020
[28,] 131.3049360 -278.2387397
[29,] 25.3716615 131.3049360
[30,] -166.5693257 25.3716615
[31,] -234.9974879 -166.5693257
[32,] 179.1153371 -234.9974879
[33,] -4.2129132 179.1153371
[34,] 345.9973115 -4.2129132
[35,] 410.6408990 345.9973115
[36,] 251.9126487 410.6408990
[37,] 0.9691494 251.9126487
[38,] 380.5230498 0.9691494
[39,] 242.1845747 380.5230498
[40,] 305.5614367 242.1845747
[41,] 422.1537240 305.5614367
[42,] 505.8665491 422.1537240
[43,] 954.5666373 505.8665491
[44,] 368.6794623 954.5666373
[45,] 324.8256501 368.6794623
[46,] -135.6179374 324.8256501
[47,] 356.4615249 -135.6179374
[48,] -98.6129132 356.4615249
[49,] 162.8794623 -98.6129132
[50,] 145.0871750 162.8794623
[51,] 75.0948876 145.0871750
[52,] 187.2923755 75.0948876
[53,] 124.7949758 187.2923755
[54,] 72.9822390 124.7949758
[55,] -43.5741735 72.9822390
[56,] 52.7565889 -43.5741735
[57,] -69.9178492 52.7565889
[58,] -118.0923755 -69.9178492
[59,] -305.6667255 -118.0923755
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 116.1619658 -205.1124723
2 -94.0661963 116.1619658
3 83.9415163 -94.0661963
4 -286.6816216 83.9415163
5 -271.3458350 -286.6816216
6 -365.5818862 -271.3458350
7 -540.9974879 -365.5818862
8 -529.1411636 -540.9974879
9 -467.3411636 -529.1411636
10 -137.7461878 -467.3411636
11 -392.9743499 -137.7461878
12 -357.8308506 -392.9743499
13 -84.5052888 -357.8308506
14 -348.4258265 -84.5052888
15 -122.9822390 -348.4258265
16 -337.4771266 -122.9822390
17 -300.9745263 -337.4771266
18 -46.6975761 -300.9745263
19 -134.9974879 -46.6975761
20 -71.4102248 -134.9974879
21 216.6462760 -71.4102248
22 45.4591892 216.6462760
23 -68.4613485 45.4591892
24 409.6435875 -68.4613485
25 -195.5052888 409.6435875
26 -83.1182020 -195.5052888
27 -278.2387397 -83.1182020
28 131.3049360 -278.2387397
29 25.3716615 131.3049360
30 -166.5693257 25.3716615
31 -234.9974879 -166.5693257
32 179.1153371 -234.9974879
33 -4.2129132 179.1153371
34 345.9973115 -4.2129132
35 410.6408990 345.9973115
36 251.9126487 410.6408990
37 0.9691494 251.9126487
38 380.5230498 0.9691494
39 242.1845747 380.5230498
40 305.5614367 242.1845747
41 422.1537240 305.5614367
42 505.8665491 422.1537240
43 954.5666373 505.8665491
44 368.6794623 954.5666373
45 324.8256501 368.6794623
46 -135.6179374 324.8256501
47 356.4615249 -135.6179374
48 -98.6129132 356.4615249
49 162.8794623 -98.6129132
50 145.0871750 162.8794623
51 75.0948876 145.0871750
52 187.2923755 75.0948876
53 124.7949758 187.2923755
54 72.9822390 124.7949758
55 -43.5741735 72.9822390
56 52.7565889 -43.5741735
57 -69.9178492 52.7565889
58 -118.0923755 -69.9178492
59 -305.6667255 -118.0923755
> 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/7srhb1258735161.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/8bqdx1258735161.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/9csq31258735161.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/10tglf1258735161.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/11quwb1258735161.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/12xq961258735161.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/13zw7e1258735161.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/14pg3g1258735161.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/15ayvp1258735161.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/166rv31258735161.tab")
+ }
> system("convert tmp/1hrl11258735161.ps tmp/1hrl11258735161.png")
> system("convert tmp/21xzq1258735161.ps tmp/21xzq1258735161.png")
> system("convert tmp/3a3mq1258735161.ps tmp/3a3mq1258735161.png")
> system("convert tmp/47u061258735161.ps tmp/47u061258735161.png")
> system("convert tmp/58qcr1258735161.ps tmp/58qcr1258735161.png")
> system("convert tmp/6gki61258735161.ps tmp/6gki61258735161.png")
> system("convert tmp/7srhb1258735161.ps tmp/7srhb1258735161.png")
> system("convert tmp/8bqdx1258735161.ps tmp/8bqdx1258735161.png")
> system("convert tmp/9csq31258735161.ps tmp/9csq31258735161.png")
> system("convert tmp/10tglf1258735161.ps tmp/10tglf1258735161.png")
>
>
> proc.time()
user system elapsed
2.392 1.564 3.425