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(6802.96,0,7132.68,0,7073.29,0,7264.5,0,7105.33,0,7218.71,0,7225.72,0,7354.25,0,7745.46,0,8070.26,0,8366.33,0,8667.51,0,8854.34,0,9218.1,0,9332.9,0,9358.31,0,9248.66,0,9401.2,0,9652.04,0,9957.38,0,10110.63,0,10169.26,0,10343.78,0,10750.21,0,11337.5,0,11786.96,0,12083.04,0,12007.74,0,11745.93,0,11051.51,0,11445.9,0,11924.88,0,12247.63,0,12690.91,0,12910.7,0,13202.12,0,13654.67,0,13862.82,0,13523.93,0,14211.17,0,14510.35,0,14289.23,0,14111.82,0,13086.59,0,13351.54,0,13747.69,0,12855.61,0,12926.93,0,12121.95,1,11731.65,1,11639.51,1,12163.78,1,12029.53,1,11234.18,1,9852.13,1,9709.04,1,9332.75,1,7108.6,1,6691.49,1,6143.05,1),dim=c(2,60),dimnames=list(c('Y','X'),1:60))
> y <- array(NA,dim=c(2,60),dimnames=list(c('Y','X'),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 = '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
Y X M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 6802.96 0 1 0 0 0 0 0 0 0 0 0 0 1
2 7132.68 0 0 1 0 0 0 0 0 0 0 0 0 2
3 7073.29 0 0 0 1 0 0 0 0 0 0 0 0 3
4 7264.50 0 0 0 0 1 0 0 0 0 0 0 0 4
5 7105.33 0 0 0 0 0 1 0 0 0 0 0 0 5
6 7218.71 0 0 0 0 0 0 1 0 0 0 0 0 6
7 7225.72 0 0 0 0 0 0 0 1 0 0 0 0 7
8 7354.25 0 0 0 0 0 0 0 0 1 0 0 0 8
9 7745.46 0 0 0 0 0 0 0 0 0 1 0 0 9
10 8070.26 0 0 0 0 0 0 0 0 0 0 1 0 10
11 8366.33 0 0 0 0 0 0 0 0 0 0 0 1 11
12 8667.51 0 0 0 0 0 0 0 0 0 0 0 0 12
13 8854.34 0 1 0 0 0 0 0 0 0 0 0 0 13
14 9218.10 0 0 1 0 0 0 0 0 0 0 0 0 14
15 9332.90 0 0 0 1 0 0 0 0 0 0 0 0 15
16 9358.31 0 0 0 0 1 0 0 0 0 0 0 0 16
17 9248.66 0 0 0 0 0 1 0 0 0 0 0 0 17
18 9401.20 0 0 0 0 0 0 1 0 0 0 0 0 18
19 9652.04 0 0 0 0 0 0 0 1 0 0 0 0 19
20 9957.38 0 0 0 0 0 0 0 0 1 0 0 0 20
21 10110.63 0 0 0 0 0 0 0 0 0 1 0 0 21
22 10169.26 0 0 0 0 0 0 0 0 0 0 1 0 22
23 10343.78 0 0 0 0 0 0 0 0 0 0 0 1 23
24 10750.21 0 0 0 0 0 0 0 0 0 0 0 0 24
25 11337.50 0 1 0 0 0 0 0 0 0 0 0 0 25
26 11786.96 0 0 1 0 0 0 0 0 0 0 0 0 26
27 12083.04 0 0 0 1 0 0 0 0 0 0 0 0 27
28 12007.74 0 0 0 0 1 0 0 0 0 0 0 0 28
29 11745.93 0 0 0 0 0 1 0 0 0 0 0 0 29
30 11051.51 0 0 0 0 0 0 1 0 0 0 0 0 30
31 11445.90 0 0 0 0 0 0 0 1 0 0 0 0 31
32 11924.88 0 0 0 0 0 0 0 0 1 0 0 0 32
33 12247.63 0 0 0 0 0 0 0 0 0 1 0 0 33
34 12690.91 0 0 0 0 0 0 0 0 0 0 1 0 34
35 12910.70 0 0 0 0 0 0 0 0 0 0 0 1 35
36 13202.12 0 0 0 0 0 0 0 0 0 0 0 0 36
37 13654.67 0 1 0 0 0 0 0 0 0 0 0 0 37
38 13862.82 0 0 1 0 0 0 0 0 0 0 0 0 38
39 13523.93 0 0 0 1 0 0 0 0 0 0 0 0 39
40 14211.17 0 0 0 0 1 0 0 0 0 0 0 0 40
41 14510.35 0 0 0 0 0 1 0 0 0 0 0 0 41
42 14289.23 0 0 0 0 0 0 1 0 0 0 0 0 42
43 14111.82 0 0 0 0 0 0 0 1 0 0 0 0 43
44 13086.59 0 0 0 0 0 0 0 0 1 0 0 0 44
45 13351.54 0 0 0 0 0 0 0 0 0 1 0 0 45
46 13747.69 0 0 0 0 0 0 0 0 0 0 1 0 46
47 12855.61 0 0 0 0 0 0 0 0 0 0 0 1 47
48 12926.93 0 0 0 0 0 0 0 0 0 0 0 0 48
49 12121.95 1 1 0 0 0 0 0 0 0 0 0 0 49
50 11731.65 1 0 1 0 0 0 0 0 0 0 0 0 50
51 11639.51 1 0 0 1 0 0 0 0 0 0 0 0 51
52 12163.78 1 0 0 0 1 0 0 0 0 0 0 0 52
53 12029.53 1 0 0 0 0 1 0 0 0 0 0 0 53
54 11234.18 1 0 0 0 0 0 1 0 0 0 0 0 54
55 9852.13 1 0 0 0 0 0 0 1 0 0 0 0 55
56 9709.04 1 0 0 0 0 0 0 0 1 0 0 0 56
57 9332.75 1 0 0 0 0 0 0 0 0 1 0 0 57
58 7108.60 1 0 0 0 0 0 0 0 0 0 1 0 58
59 6691.49 1 0 0 0 0 0 0 0 0 0 0 1 59
60 6143.05 1 0 0 0 0 0 0 0 0 0 0 0 60
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X M1 M2 M3 M4
5272.09 -5964.88 2128.75 2147.05 1957.28 2053.99
M5 M6 M7 M8 M9 M10
1806.99 1344.14 988.84 763.89 741.21 367.09
M11 t
69.47 173.86
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-3595.6 -627.7 -125.3 545.9 2167.0
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5272.09 677.48 7.782 6.21e-10 ***
X -5964.88 556.88 -10.711 4.38e-14 ***
M1 2128.75 785.03 2.712 0.00938 **
M2 2147.05 782.72 2.743 0.00865 **
M3 1957.28 780.63 2.507 0.01576 *
M4 2053.99 778.75 2.638 0.01135 *
M5 1806.99 777.09 2.325 0.02452 *
M6 1344.14 775.65 1.733 0.08981 .
M7 988.84 774.42 1.277 0.20805
M8 763.89 773.42 0.988 0.32848
M9 741.21 772.64 0.959 0.34241
M10 367.09 772.08 0.475 0.63671
M11 69.47 771.75 0.090 0.92866
t 173.86 13.13 13.245 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1220 on 46 degrees of freedom
Multiple R-squared: 0.7975, Adjusted R-squared: 0.7403
F-statistic: 13.94 on 13 and 46 DF, p-value: 7.491e-12
> 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,] 1.797658e-04 3.595315e-04 0.9998202
[2,] 9.202223e-06 1.840445e-05 0.9999908
[3,] 8.107991e-06 1.621598e-05 0.9999919
[4,] 8.694427e-06 1.738885e-05 0.9999913
[5,] 1.015456e-06 2.030913e-06 0.9999990
[6,] 1.183152e-07 2.366304e-07 0.9999999
[7,] 2.318936e-08 4.637872e-08 1.0000000
[8,] 2.820750e-09 5.641499e-09 1.0000000
[9,] 6.075865e-10 1.215173e-09 1.0000000
[10,] 1.845319e-10 3.690638e-10 1.0000000
[11,] 3.215567e-10 6.431135e-10 1.0000000
[12,] 9.075300e-11 1.815060e-10 1.0000000
[13,] 2.405824e-11 4.811647e-11 1.0000000
[14,] 1.150143e-09 2.300287e-09 1.0000000
[15,] 1.510493e-09 3.020986e-09 1.0000000
[16,] 5.408763e-10 1.081753e-09 1.0000000
[17,] 2.170964e-10 4.341927e-10 1.0000000
[18,] 5.724864e-11 1.144973e-10 1.0000000
[19,] 1.086754e-11 2.173508e-11 1.0000000
[20,] 1.625667e-12 3.251334e-12 1.0000000
[21,] 1.043869e-12 2.087738e-12 1.0000000
[22,] 4.435123e-13 8.870245e-13 1.0000000
[23,] 3.469823e-12 6.939647e-12 1.0000000
[24,] 1.078932e-11 2.157864e-11 1.0000000
[25,] 5.439909e-10 1.087982e-09 1.0000000
[26,] 9.441249e-09 1.888250e-08 1.0000000
[27,] 4.665443e-09 9.330886e-09 1.0000000
> postscript(file="/var/www/html/rcomp/tmp/1iog21258740313.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/2nseb1258740313.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/3maj61258740313.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/4wol91258740313.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/5rz7n1258740313.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
-771.73300 -634.17100 -677.65300 -757.00900 -843.03900 -440.66500
7 8 9 10 11 12
-252.21100 -72.58700 167.44900 692.50700 1112.33900 1309.13700
13 14 15 16 17 18
-806.63683 -635.03483 -504.32683 -749.48283 -785.99283 -344.45883
19 20 21 22 23 24
87.82517 444.25917 446.33517 705.22317 1003.50517 1305.55317
25 26 27 28 29 30
-409.76067 -152.45867 159.52933 -186.33667 -375.00667 -780.43267
31 32 33 34 35 36
-204.59867 325.47533 497.05133 1140.58933 1484.14133 1671.17933
37 38 39 40 41 42
-178.87450 -162.88250 -485.86450 -69.19050 303.12950 371.00350
43 44 45 46 47 48
375.03750 -599.09850 -485.32250 111.08550 -657.23250 -690.29450
49 50 51 52 53 54
2167.00500 1584.54700 1508.31500 1762.01900 1700.90900 1194.55300
55 56 57 58 59 60
-6.05300 -98.04900 -625.51300 -2649.40500 -2942.75300 -3595.57500
> postscript(file="/var/www/html/rcomp/tmp/6iamd1258740313.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 -771.73300 NA
1 -634.17100 -771.73300
2 -677.65300 -634.17100
3 -757.00900 -677.65300
4 -843.03900 -757.00900
5 -440.66500 -843.03900
6 -252.21100 -440.66500
7 -72.58700 -252.21100
8 167.44900 -72.58700
9 692.50700 167.44900
10 1112.33900 692.50700
11 1309.13700 1112.33900
12 -806.63683 1309.13700
13 -635.03483 -806.63683
14 -504.32683 -635.03483
15 -749.48283 -504.32683
16 -785.99283 -749.48283
17 -344.45883 -785.99283
18 87.82517 -344.45883
19 444.25917 87.82517
20 446.33517 444.25917
21 705.22317 446.33517
22 1003.50517 705.22317
23 1305.55317 1003.50517
24 -409.76067 1305.55317
25 -152.45867 -409.76067
26 159.52933 -152.45867
27 -186.33667 159.52933
28 -375.00667 -186.33667
29 -780.43267 -375.00667
30 -204.59867 -780.43267
31 325.47533 -204.59867
32 497.05133 325.47533
33 1140.58933 497.05133
34 1484.14133 1140.58933
35 1671.17933 1484.14133
36 -178.87450 1671.17933
37 -162.88250 -178.87450
38 -485.86450 -162.88250
39 -69.19050 -485.86450
40 303.12950 -69.19050
41 371.00350 303.12950
42 375.03750 371.00350
43 -599.09850 375.03750
44 -485.32250 -599.09850
45 111.08550 -485.32250
46 -657.23250 111.08550
47 -690.29450 -657.23250
48 2167.00500 -690.29450
49 1584.54700 2167.00500
50 1508.31500 1584.54700
51 1762.01900 1508.31500
52 1700.90900 1762.01900
53 1194.55300 1700.90900
54 -6.05300 1194.55300
55 -98.04900 -6.05300
56 -625.51300 -98.04900
57 -2649.40500 -625.51300
58 -2942.75300 -2649.40500
59 -3595.57500 -2942.75300
60 NA -3595.57500
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -634.17100 -771.73300
[2,] -677.65300 -634.17100
[3,] -757.00900 -677.65300
[4,] -843.03900 -757.00900
[5,] -440.66500 -843.03900
[6,] -252.21100 -440.66500
[7,] -72.58700 -252.21100
[8,] 167.44900 -72.58700
[9,] 692.50700 167.44900
[10,] 1112.33900 692.50700
[11,] 1309.13700 1112.33900
[12,] -806.63683 1309.13700
[13,] -635.03483 -806.63683
[14,] -504.32683 -635.03483
[15,] -749.48283 -504.32683
[16,] -785.99283 -749.48283
[17,] -344.45883 -785.99283
[18,] 87.82517 -344.45883
[19,] 444.25917 87.82517
[20,] 446.33517 444.25917
[21,] 705.22317 446.33517
[22,] 1003.50517 705.22317
[23,] 1305.55317 1003.50517
[24,] -409.76067 1305.55317
[25,] -152.45867 -409.76067
[26,] 159.52933 -152.45867
[27,] -186.33667 159.52933
[28,] -375.00667 -186.33667
[29,] -780.43267 -375.00667
[30,] -204.59867 -780.43267
[31,] 325.47533 -204.59867
[32,] 497.05133 325.47533
[33,] 1140.58933 497.05133
[34,] 1484.14133 1140.58933
[35,] 1671.17933 1484.14133
[36,] -178.87450 1671.17933
[37,] -162.88250 -178.87450
[38,] -485.86450 -162.88250
[39,] -69.19050 -485.86450
[40,] 303.12950 -69.19050
[41,] 371.00350 303.12950
[42,] 375.03750 371.00350
[43,] -599.09850 375.03750
[44,] -485.32250 -599.09850
[45,] 111.08550 -485.32250
[46,] -657.23250 111.08550
[47,] -690.29450 -657.23250
[48,] 2167.00500 -690.29450
[49,] 1584.54700 2167.00500
[50,] 1508.31500 1584.54700
[51,] 1762.01900 1508.31500
[52,] 1700.90900 1762.01900
[53,] 1194.55300 1700.90900
[54,] -6.05300 1194.55300
[55,] -98.04900 -6.05300
[56,] -625.51300 -98.04900
[57,] -2649.40500 -625.51300
[58,] -2942.75300 -2649.40500
[59,] -3595.57500 -2942.75300
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -634.17100 -771.73300
2 -677.65300 -634.17100
3 -757.00900 -677.65300
4 -843.03900 -757.00900
5 -440.66500 -843.03900
6 -252.21100 -440.66500
7 -72.58700 -252.21100
8 167.44900 -72.58700
9 692.50700 167.44900
10 1112.33900 692.50700
11 1309.13700 1112.33900
12 -806.63683 1309.13700
13 -635.03483 -806.63683
14 -504.32683 -635.03483
15 -749.48283 -504.32683
16 -785.99283 -749.48283
17 -344.45883 -785.99283
18 87.82517 -344.45883
19 444.25917 87.82517
20 446.33517 444.25917
21 705.22317 446.33517
22 1003.50517 705.22317
23 1305.55317 1003.50517
24 -409.76067 1305.55317
25 -152.45867 -409.76067
26 159.52933 -152.45867
27 -186.33667 159.52933
28 -375.00667 -186.33667
29 -780.43267 -375.00667
30 -204.59867 -780.43267
31 325.47533 -204.59867
32 497.05133 325.47533
33 1140.58933 497.05133
34 1484.14133 1140.58933
35 1671.17933 1484.14133
36 -178.87450 1671.17933
37 -162.88250 -178.87450
38 -485.86450 -162.88250
39 -69.19050 -485.86450
40 303.12950 -69.19050
41 371.00350 303.12950
42 375.03750 371.00350
43 -599.09850 375.03750
44 -485.32250 -599.09850
45 111.08550 -485.32250
46 -657.23250 111.08550
47 -690.29450 -657.23250
48 2167.00500 -690.29450
49 1584.54700 2167.00500
50 1508.31500 1584.54700
51 1762.01900 1508.31500
52 1700.90900 1762.01900
53 1194.55300 1700.90900
54 -6.05300 1194.55300
55 -98.04900 -6.05300
56 -625.51300 -98.04900
57 -2649.40500 -625.51300
58 -2942.75300 -2649.40500
59 -3595.57500 -2942.75300
> 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/72r8p1258740313.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/85a621258740313.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/98zz41258740313.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/103na11258740313.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/11zqgp1258740313.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/12dcj01258740313.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/138vaf1258740313.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/14qekf1258740313.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/15p9jq1258740313.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/16wo001258740313.tab")
+ }
> system("convert tmp/1iog21258740313.ps tmp/1iog21258740313.png")
> system("convert tmp/2nseb1258740313.ps tmp/2nseb1258740313.png")
> system("convert tmp/3maj61258740313.ps tmp/3maj61258740313.png")
> system("convert tmp/4wol91258740313.ps tmp/4wol91258740313.png")
> system("convert tmp/5rz7n1258740313.ps tmp/5rz7n1258740313.png")
> system("convert tmp/6iamd1258740313.ps tmp/6iamd1258740313.png")
> system("convert tmp/72r8p1258740313.ps tmp/72r8p1258740313.png")
> system("convert tmp/85a621258740313.ps tmp/85a621258740313.png")
> system("convert tmp/98zz41258740313.ps tmp/98zz41258740313.png")
> system("convert tmp/103na11258740313.ps tmp/103na11258740313.png")
>
>
> proc.time()
user system elapsed
2.385 1.545 2.791