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(-999.00
+ ,-999.00
+ ,38.60
+ ,6654.00
+ ,5712.00
+ ,645.00
+ ,3.00
+ ,5.00
+ ,3.00
+ ,3.30
+ ,6.30
+ ,2.00
+ ,4.50
+ ,1.00
+ ,6600.00
+ ,42.00
+ ,3.00
+ ,1.00
+ ,3.00
+ ,8.30
+ ,-999.00
+ ,-999.00
+ ,14.00
+ ,3.39
+ ,44.50
+ ,60.00
+ ,1.00
+ ,1.00
+ ,1.00
+ ,12.50
+ ,-999.00
+ ,-999.00
+ ,-999.00
+ ,0.92
+ ,5.70
+ ,25.00
+ ,5.00
+ ,2.00
+ ,3.00
+ ,16.50
+ ,2.10
+ ,1.80
+ ,69.00
+ ,2547.00
+ ,4603.00
+ ,624.00
+ ,3.00
+ ,5.00
+ ,4.00
+ ,3.90
+ ,9.10
+ ,0.70
+ ,27.00
+ ,10.55
+ ,179.50
+ ,180.00
+ ,4.00
+ ,4.00
+ ,4.00
+ ,9.80
+ ,15.80
+ ,3.90
+ ,19.00
+ ,0.02
+ ,0.30
+ ,35.00
+ ,1.00
+ ,1.00
+ ,1.00
+ ,19.70
+ ,5.20
+ ,1.00
+ ,30.40
+ ,160.00
+ ,169.00
+ ,392.00
+ ,4.00
+ ,5.00
+ ,4.00
+ ,6.20
+ ,10.90
+ ,3.60
+ ,28.00
+ ,3.30
+ ,25.60
+ ,63.00
+ ,1.00
+ ,2.00
+ ,1.00
+ ,14.50
+ ,8.30
+ ,1.40
+ ,50.00
+ ,52.16
+ ,440.00
+ ,230.00
+ ,1.00
+ ,1.00
+ ,1.00
+ ,9.70
+ ,11.00
+ ,1.50
+ ,7.00
+ ,0.43
+ ,6.40
+ ,112.00
+ ,5.00
+ ,4.00
+ ,4.00
+ ,12.50
+ ,3.20
+ ,0.70
+ ,30.00
+ ,465.00
+ ,423.00
+ ,281.00
+ ,5.00
+ ,5.00
+ ,5.00
+ ,3.90
+ ,7.60
+ ,2.70
+ ,-999.00
+ ,0.55
+ ,2.40
+ ,-999.00
+ ,2.00
+ ,1.00
+ ,2.00
+ ,10.30
+ ,-999.00
+ ,-999.00
+ ,40.00
+ ,187.10
+ ,419.00
+ ,365.00
+ ,5.00
+ ,5.00
+ ,5.00
+ ,3.10
+ ,6.30
+ ,2.10
+ ,3.50
+ ,0.08
+ ,1.20
+ ,42.00
+ ,1.00
+ ,1.00
+ ,1.00
+ ,8.40
+ ,8.60
+ ,0.00
+ ,50.00
+ ,3.00
+ ,25.00
+ ,28.00
+ ,2.00
+ ,2.00
+ ,2.00
+ ,8.60
+ ,6.60
+ ,4.10
+ ,6.00
+ ,0.79
+ ,3500.00
+ ,42.00
+ ,2.00
+ ,2.00
+ ,2.00
+ ,10.70
+ ,9.50
+ ,1.20
+ ,10.40
+ ,0.20
+ ,5.00
+ ,120.00
+ ,2.00
+ ,2.00
+ ,2.00
+ ,10.70
+ ,4.80
+ ,1.30
+ ,34.00
+ ,1.41
+ ,17.50
+ ,-999.00
+ ,1.00
+ ,2.00
+ ,1.00
+ ,6.10
+ ,12.00
+ ,6.10
+ ,7.00
+ ,60.00
+ ,81.00
+ ,-999.00
+ ,1.00
+ ,1.00
+ ,1.00
+ ,18.10
+ ,-999.00
+ ,0.30
+ ,28.00
+ ,529.00
+ ,680.00
+ ,400.00
+ ,5.00
+ ,5.00
+ ,5.00
+ ,-999.00
+ ,3.30
+ ,0.50
+ ,20.00
+ ,27.66
+ ,115.00
+ ,148.00
+ ,5.00
+ ,5.00
+ ,5.00
+ ,3.80
+ ,11.00
+ ,3.40
+ ,3.90
+ ,0.12
+ ,1.00
+ ,16.00
+ ,3.00
+ ,1.00
+ ,2.00
+ ,14.40
+ ,-999.00
+ ,-999.00
+ ,39.30
+ ,207.00
+ ,406.00
+ ,252.00
+ ,1.00
+ ,4.00
+ ,1.00
+ ,12.00
+ ,4.70
+ ,1.50
+ ,41.00
+ ,85.00
+ ,325.00
+ ,310.00
+ ,1.00
+ ,3.00
+ ,1.00
+ ,6.20
+ ,-999.00
+ ,-999.00
+ ,16.20
+ ,36.33
+ ,119.50
+ ,63.00
+ ,1.00
+ ,1.00
+ ,1.00
+ ,13.00
+ ,10.40
+ ,3.40
+ ,9.00
+ ,0.10
+ ,4.00
+ ,28.00
+ ,5.00
+ ,1.00
+ ,3.00
+ ,13.80
+ ,7.40
+ ,0.80
+ ,7.60
+ ,1.04
+ ,5.50
+ ,68.00
+ ,5.00
+ ,3.00
+ ,4.00
+ ,8.20
+ ,2.10
+ ,0.80
+ ,46.00
+ ,521.00
+ ,655.00
+ ,336.00
+ ,5.00
+ ,5.00
+ ,5.00
+ ,2.90
+ ,-999.00
+ ,-999.00
+ ,22.40
+ ,100.00
+ ,157.00
+ ,100.00
+ ,1.00
+ ,1.00
+ ,1.00
+ ,10.80
+ ,-999.00
+ ,-999.00
+ ,16.30
+ ,35.00
+ ,56.00
+ ,33.00
+ ,3.00
+ ,5.00
+ ,4.00
+ ,-999.00
+ ,7.70
+ ,1.40
+ ,2.60
+ ,0.01
+ ,0.14
+ ,21.50
+ ,5.00
+ ,2.00
+ ,4.00
+ ,9.10
+ ,17.90
+ ,2.00
+ ,24.00
+ ,0.01
+ ,0.25
+ ,50.00
+ ,1.00
+ ,1.00
+ ,1.00
+ ,19.90
+ ,6.10
+ ,1.90
+ ,100.00
+ ,62.00
+ ,1320.00
+ ,267.00
+ ,1.00
+ ,1.00
+ ,1.00
+ ,8.00
+ ,8.20
+ ,2.40
+ ,-999.00
+ ,0.12
+ ,3.00
+ ,30.00
+ ,2.00
+ ,1.00
+ ,1.00
+ ,10.60
+ ,8.40
+ ,2.80
+ ,-999.00
+ ,1.35
+ ,8.10
+ ,45.00
+ ,3.00
+ ,1.00
+ ,3.00
+ ,11.20
+ ,11.90
+ ,1.30
+ ,3.20
+ ,0.02
+ ,0.40
+ ,19.00
+ ,4.00
+ ,1.00
+ ,3.00
+ ,13.20
+ ,10.80
+ ,2.00
+ ,2.00
+ ,0.05
+ ,0.33
+ ,30.00
+ ,4.00
+ ,1.00
+ ,3.00
+ ,12.80
+ ,13.80
+ ,5.60
+ ,5.00
+ ,1.70
+ ,6.30
+ ,12.00
+ ,2.00
+ ,1.00
+ ,1.00
+ ,19.40
+ ,14.30
+ ,3.10
+ ,6.50
+ ,3.50
+ ,10.80
+ ,120.00
+ ,2.00
+ ,1.00
+ ,1.00
+ ,17.40
+ ,-999.00
+ ,1.00
+ ,23.60
+ ,250.00
+ ,490.00
+ ,440.00
+ ,5.00
+ ,5.00
+ ,5.00
+ ,-999.00
+ ,15.20
+ ,1.80
+ ,12.00
+ ,0.48
+ ,15.50
+ ,140.00
+ ,2.00
+ ,2.00
+ ,2.00
+ ,17.00
+ ,10.00
+ ,0.90
+ ,20.20
+ ,10.00
+ ,115.00
+ ,170.00
+ ,4.00
+ ,4.00
+ ,4.00
+ ,10.90
+ ,11.90
+ ,1.80
+ ,13.00
+ ,1.62
+ ,11.40
+ ,17.00
+ ,2.00
+ ,1.00
+ ,2.00
+ ,13.70
+ ,6.50
+ ,1.90
+ ,27.00
+ ,192.00
+ ,180.00
+ ,115.00
+ ,4.00
+ ,4.00
+ ,4.00
+ ,8.40
+ ,7.50
+ ,0.90
+ ,18.00
+ ,2.50
+ ,12.10
+ ,31.00
+ ,5.00
+ ,5.00
+ ,5.00
+ ,8.40
+ ,-999.00
+ ,-999.00
+ ,13.70
+ ,4.29
+ ,39.20
+ ,63.00
+ ,2.00
+ ,2.00
+ ,2.00
+ ,12.50
+ ,10.60
+ ,2.60
+ ,4.70
+ ,0.28
+ ,1.90
+ ,21.00
+ ,3.00
+ ,1.00
+ ,3.00
+ ,13.20
+ ,7.40
+ ,2.40
+ ,9.80
+ ,4.24
+ ,50.40
+ ,52.00
+ ,1.00
+ ,1.00
+ ,1.00
+ ,9.80
+ ,8.40
+ ,1.20
+ ,29.00
+ ,6.80
+ ,179.00
+ ,164.00
+ ,2.00
+ ,3.00
+ ,2.00
+ ,9.60
+ ,5.70
+ ,0.90
+ ,7.00
+ ,0.75
+ ,12.30
+ ,225.00
+ ,2.00
+ ,2.00
+ ,2.00
+ ,6.60
+ ,4.90
+ ,0.50
+ ,6.00
+ ,3.60
+ ,21.00
+ ,225.00
+ ,3.00
+ ,2.00
+ ,3.00
+ ,5.40
+ ,-999.00
+ ,-999.00
+ ,17.00
+ ,14.83
+ ,98.20
+ ,150.00
+ ,5.00
+ ,5.00
+ ,5.00
+ ,2.60
+ ,3.20
+ ,0.60
+ ,20.00
+ ,55.50
+ ,175.00
+ ,151.00
+ ,5.00
+ ,5.00
+ ,5.00
+ ,3.80
+ ,-999.00
+ ,-999.00
+ ,12.70
+ ,1.40
+ ,12.50
+ ,90.00
+ ,2.00
+ ,2.00
+ ,2.00
+ ,11.00
+ ,8.10
+ ,2.20
+ ,3.50
+ ,0.06
+ ,1.00
+ ,-999.00
+ ,3.00
+ ,1.00
+ ,2.00
+ ,10.30
+ ,11.00
+ ,2.30
+ ,4.50
+ ,0.90
+ ,2.60
+ ,60.00
+ ,2.00
+ ,1.00
+ ,2.00
+ ,13.30
+ ,4.90
+ ,0.50
+ ,7.50
+ ,2.00
+ ,12.30
+ ,200.00
+ ,3.00
+ ,1.00
+ ,3.00
+ ,5.40
+ ,13.20
+ ,2.60
+ ,2.30
+ ,0.10
+ ,2.50
+ ,46.00
+ ,3.00
+ ,2.00
+ ,2.00
+ ,15.80
+ ,9.70
+ ,0.60
+ ,24.00
+ ,4.19
+ ,58.00
+ ,210.00
+ ,4.00
+ ,3.00
+ ,4.00
+ ,10.30
+ ,12.80
+ ,6.60
+ ,3.00
+ ,3.50
+ ,3.90
+ ,14.00
+ ,2.00
+ ,1.00
+ ,1.00
+ ,19.40
+ ,-999.00
+ ,-999.00
+ ,13.00
+ ,4.05
+ ,17.00
+ ,38.00
+ ,3.00
+ ,1.00
+ ,1.00
+ ,-999.00)
+ ,dim=c(10
+ ,62)
+ ,dimnames=list(c('SWS'
+ ,'PS'
+ ,'L'
+ ,'Wb'
+ ,'Wbr'
+ ,'Tg'
+ ,'P'
+ ,'S'
+ ,'D'
+ ,'TS')
+ ,1:62))
> y <- array(NA,dim=c(10,62),dimnames=list(c('SWS','PS','L','Wb','Wbr','Tg','P','S','D','TS'),1:62))
> for (i in 1:dim(x)[1])
+ {
+ for (j in 1:dim(x)[2])
+ {
+ y[i,j] <- as.numeric(x[i,j])
+ }
+ }
> par20 = ''
> par19 = ''
> par18 = ''
> par17 = ''
> par16 = ''
> par15 = ''
> par14 = ''
> par13 = ''
> par12 = ''
> par11 = ''
> par10 = ''
> par9 = ''
> par8 = ''
> par7 = ''
> par6 = ''
> par5 = ''
> par4 = ''
> par3 = 'No Linear Trend'
> par2 = 'Do not include Seasonal Dummies'
> par1 = '10'
> ylab = ''
> xlab = ''
> main = ''
> #'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
TS SWS PS L Wb Wbr Tg P S D
1 3.3 -999.0 -999.0 38.6 6654.00 5712.00 645.0 3 5 3
2 8.3 6.3 2.0 4.5 1.00 6600.00 42.0 3 1 3
3 12.5 -999.0 -999.0 14.0 3.39 44.50 60.0 1 1 1
4 16.5 -999.0 -999.0 -999.0 0.92 5.70 25.0 5 2 3
5 3.9 2.1 1.8 69.0 2547.00 4603.00 624.0 3 5 4
6 9.8 9.1 0.7 27.0 10.55 179.50 180.0 4 4 4
7 19.7 15.8 3.9 19.0 0.02 0.30 35.0 1 1 1
8 6.2 5.2 1.0 30.4 160.00 169.00 392.0 4 5 4
9 14.5 10.9 3.6 28.0 3.30 25.60 63.0 1 2 1
10 9.7 8.3 1.4 50.0 52.16 440.00 230.0 1 1 1
11 12.5 11.0 1.5 7.0 0.43 6.40 112.0 5 4 4
12 3.9 3.2 0.7 30.0 465.00 423.00 281.0 5 5 5
13 10.3 7.6 2.7 -999.0 0.55 2.40 -999.0 2 1 2
14 3.1 -999.0 -999.0 40.0 187.10 419.00 365.0 5 5 5
15 8.4 6.3 2.1 3.5 0.08 1.20 42.0 1 1 1
16 8.6 8.6 0.0 50.0 3.00 25.00 28.0 2 2 2
17 10.7 6.6 4.1 6.0 0.79 3500.00 42.0 2 2 2
18 10.7 9.5 1.2 10.4 0.20 5.00 120.0 2 2 2
19 6.1 4.8 1.3 34.0 1.41 17.50 -999.0 1 2 1
20 18.1 12.0 6.1 7.0 60.00 81.00 -999.0 1 1 1
21 -999.0 -999.0 0.3 28.0 529.00 680.00 400.0 5 5 5
22 3.8 3.3 0.5 20.0 27.66 115.00 148.0 5 5 5
23 14.4 11.0 3.4 3.9 0.12 1.00 16.0 3 1 2
24 12.0 -999.0 -999.0 39.3 207.00 406.00 252.0 1 4 1
25 6.2 4.7 1.5 41.0 85.00 325.00 310.0 1 3 1
26 13.0 -999.0 -999.0 16.2 36.33 119.50 63.0 1 1 1
27 13.8 10.4 3.4 9.0 0.10 4.00 28.0 5 1 3
28 8.2 7.4 0.8 7.6 1.04 5.50 68.0 5 3 4
29 2.9 2.1 0.8 46.0 521.00 655.00 336.0 5 5 5
30 10.8 -999.0 -999.0 22.4 100.00 157.00 100.0 1 1 1
31 -999.0 -999.0 -999.0 16.3 35.00 56.00 33.0 3 5 4
32 9.1 7.7 1.4 2.6 0.01 0.14 21.5 5 2 4
33 19.9 17.9 2.0 24.0 0.01 0.25 50.0 1 1 1
34 8.0 6.1 1.9 100.0 62.00 1320.00 267.0 1 1 1
35 10.6 8.2 2.4 -999.0 0.12 3.00 30.0 2 1 1
36 11.2 8.4 2.8 -999.0 1.35 8.10 45.0 3 1 3
37 13.2 11.9 1.3 3.2 0.02 0.40 19.0 4 1 3
38 12.8 10.8 2.0 2.0 0.05 0.33 30.0 4 1 3
39 19.4 13.8 5.6 5.0 1.70 6.30 12.0 2 1 1
40 17.4 14.3 3.1 6.5 3.50 10.80 120.0 2 1 1
41 -999.0 -999.0 1.0 23.6 250.00 490.00 440.0 5 5 5
42 17.0 15.2 1.8 12.0 0.48 15.50 140.0 2 2 2
43 10.9 10.0 0.9 20.2 10.00 115.00 170.0 4 4 4
44 13.7 11.9 1.8 13.0 1.62 11.40 17.0 2 1 2
45 8.4 6.5 1.9 27.0 192.00 180.00 115.0 4 4 4
46 8.4 7.5 0.9 18.0 2.50 12.10 31.0 5 5 5
47 12.5 -999.0 -999.0 13.7 4.29 39.20 63.0 2 2 2
48 13.2 10.6 2.6 4.7 0.28 1.90 21.0 3 1 3
49 9.8 7.4 2.4 9.8 4.24 50.40 52.0 1 1 1
50 9.6 8.4 1.2 29.0 6.80 179.00 164.0 2 3 2
51 6.6 5.7 0.9 7.0 0.75 12.30 225.0 2 2 2
52 5.4 4.9 0.5 6.0 3.60 21.00 225.0 3 2 3
53 2.6 -999.0 -999.0 17.0 14.83 98.20 150.0 5 5 5
54 3.8 3.2 0.6 20.0 55.50 175.00 151.0 5 5 5
55 11.0 -999.0 -999.0 12.7 1.40 12.50 90.0 2 2 2
56 10.3 8.1 2.2 3.5 0.06 1.00 -999.0 3 1 2
57 13.3 11.0 2.3 4.5 0.90 2.60 60.0 2 1 2
58 5.4 4.9 0.5 7.5 2.00 12.30 200.0 3 1 3
59 15.8 13.2 2.6 2.3 0.10 2.50 46.0 3 2 2
60 10.3 9.7 0.6 24.0 4.19 58.00 210.0 4 3 4
61 19.4 12.8 6.6 3.0 3.50 3.90 14.0 2 1 1
62 -999.0 -999.0 -999.0 13.0 4.05 17.00 38.0 3 1 1
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) SWS PS L Wb Wbr
34.36234 0.99734 -0.82684 -0.06033 0.03181 -0.00603
Tg P S D
0.05054 -36.44183 -22.24531 45.67018
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-826.62 -18.43 10.88 30.64 213.22
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 34.362342 55.213315 0.622 0.536
SWS 0.997337 0.134795 7.399 1.14e-09 ***
PS -0.826837 0.142513 -5.802 3.95e-07 ***
L -0.060330 0.096716 -0.624 0.535
Wb 0.031813 0.036817 0.864 0.392
Wbr -0.006029 0.024183 -0.249 0.804
Tg 0.050545 0.085304 0.593 0.556
P -36.441835 43.881230 -0.830 0.410
S -22.245315 29.432912 -0.756 0.453
D 45.670184 58.506465 0.781 0.439
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 176.5 on 52 degrees of freedom
Multiple R-squared: 0.5753, Adjusted R-squared: 0.5018
F-statistic: 7.826 on 9 and 52 DF, p-value: 3.506e-07
> 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,] 4.168012e-06 8.336024e-06 0.9999958
[2,] 1.255463e-07 2.510926e-07 0.9999999
[3,] 1.847229e-09 3.694459e-09 1.0000000
[4,] 4.689659e-11 9.379318e-11 1.0000000
[5,] 6.297280e-13 1.259456e-12 1.0000000
[6,] 8.103911e-15 1.620782e-14 1.0000000
[7,] 9.285437e-17 1.857087e-16 1.0000000
[8,] 8.753160e-17 1.750632e-16 1.0000000
[9,] 1.469879e-18 2.939759e-18 1.0000000
[10,] 2.275484e-20 4.550968e-20 1.0000000
[11,] 3.735718e-22 7.471436e-22 1.0000000
[12,] 8.069796e-24 1.613959e-23 1.0000000
[13,] 1.777874e-25 3.555748e-25 1.0000000
[14,] 4.181666e-27 8.363332e-27 1.0000000
[15,] 5.558873e-29 1.111775e-28 1.0000000
[16,] 8.349985e-31 1.669997e-30 1.0000000
[17,] 1.200294e-32 2.400587e-32 1.0000000
[18,] 3.205436e-34 6.410871e-34 1.0000000
[19,] 7.589995e-01 4.820011e-01 0.2410005
[20,] 7.097993e-01 5.804013e-01 0.2902007
[21,] 6.424401e-01 7.151199e-01 0.3575599
[22,] 5.866516e-01 8.266968e-01 0.4133484
[23,] 5.559664e-01 8.880672e-01 0.4440336
[24,] 6.770873e-01 6.458254e-01 0.3229127
[25,] 6.575735e-01 6.848530e-01 0.3424265
[26,] 6.968144e-01 6.063712e-01 0.3031856
[27,] 6.895709e-01 6.208582e-01 0.3104291
[28,] 7.811627e-01 4.376746e-01 0.2188373
[29,] 7.313583e-01 5.372835e-01 0.2686417
[30,] 7.239224e-01 5.521551e-01 0.2760776
[31,] 6.613951e-01 6.772098e-01 0.3386049
[32,] 5.728948e-01 8.542103e-01 0.4271052
[33,] 4.597946e-01 9.195891e-01 0.5402054
[34,] 4.408108e-01 8.816215e-01 0.5591892
[35,] 3.838696e-01 7.677392e-01 0.6161304
[36,] 2.614019e-01 5.228037e-01 0.7385981
[37,] 1.557611e-01 3.115222e-01 0.8442389
> postscript(file="/var/www/html/freestat/rcomp/tmp/1pefu1292960071.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/2pefu1292960071.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/3i5ef1292960071.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/4i5ef1292960071.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/5i5ef1292960071.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 = 62
Frequency = 1
1 2 3 4 5 6
15.295541 1.779383 159.456986 180.628282 -73.847446 12.286081
7 8 9 10 11 12
-14.800252 19.741014 6.059886 -26.380852 51.703124 7.374241
13 14 15 16 17 18
-35.400006 184.691045 -19.399271 -6.648942 18.496160 -11.525046
19 20 21 22 23 24
55.893613 39.326334 -2.947139 25.184140 11.537432 213.216987
25 26 27 28 29 30
11.958219 159.342378 38.469194 30.404768 5.356646 153.846869
31 32 33 34 35 36
-826.621428 11.305477 -18.722151 -19.334067 -42.269430 -97.203245
37 38 39 40 41 42
-1.719179 -1.073086 25.042512 15.078294 3.074535 -11.273729
43 44 45 46 47 48
12.377644 -27.311495 11.987381 31.899100 172.243632 -35.799408
49 50 51 52 53 54
-18.809349 10.454677 -17.569011 -27.628765 197.216666 24.691026
55 56 57 58 59 60
169.249548 60.618109 -29.116856 -48.521524 30.723948 -12.367897
61 62
26.573202 -778.294498
> postscript(file="/var/www/html/freestat/rcomp/tmp/6bew01292960071.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 = 62
Frequency = 1
lag(myerror, k = 1) myerror
0 15.295541 NA
1 1.779383 15.295541
2 159.456986 1.779383
3 180.628282 159.456986
4 -73.847446 180.628282
5 12.286081 -73.847446
6 -14.800252 12.286081
7 19.741014 -14.800252
8 6.059886 19.741014
9 -26.380852 6.059886
10 51.703124 -26.380852
11 7.374241 51.703124
12 -35.400006 7.374241
13 184.691045 -35.400006
14 -19.399271 184.691045
15 -6.648942 -19.399271
16 18.496160 -6.648942
17 -11.525046 18.496160
18 55.893613 -11.525046
19 39.326334 55.893613
20 -2.947139 39.326334
21 25.184140 -2.947139
22 11.537432 25.184140
23 213.216987 11.537432
24 11.958219 213.216987
25 159.342378 11.958219
26 38.469194 159.342378
27 30.404768 38.469194
28 5.356646 30.404768
29 153.846869 5.356646
30 -826.621428 153.846869
31 11.305477 -826.621428
32 -18.722151 11.305477
33 -19.334067 -18.722151
34 -42.269430 -19.334067
35 -97.203245 -42.269430
36 -1.719179 -97.203245
37 -1.073086 -1.719179
38 25.042512 -1.073086
39 15.078294 25.042512
40 3.074535 15.078294
41 -11.273729 3.074535
42 12.377644 -11.273729
43 -27.311495 12.377644
44 11.987381 -27.311495
45 31.899100 11.987381
46 172.243632 31.899100
47 -35.799408 172.243632
48 -18.809349 -35.799408
49 10.454677 -18.809349
50 -17.569011 10.454677
51 -27.628765 -17.569011
52 197.216666 -27.628765
53 24.691026 197.216666
54 169.249548 24.691026
55 60.618109 169.249548
56 -29.116856 60.618109
57 -48.521524 -29.116856
58 30.723948 -48.521524
59 -12.367897 30.723948
60 26.573202 -12.367897
61 -778.294498 26.573202
62 NA -778.294498
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 1.779383 15.295541
[2,] 159.456986 1.779383
[3,] 180.628282 159.456986
[4,] -73.847446 180.628282
[5,] 12.286081 -73.847446
[6,] -14.800252 12.286081
[7,] 19.741014 -14.800252
[8,] 6.059886 19.741014
[9,] -26.380852 6.059886
[10,] 51.703124 -26.380852
[11,] 7.374241 51.703124
[12,] -35.400006 7.374241
[13,] 184.691045 -35.400006
[14,] -19.399271 184.691045
[15,] -6.648942 -19.399271
[16,] 18.496160 -6.648942
[17,] -11.525046 18.496160
[18,] 55.893613 -11.525046
[19,] 39.326334 55.893613
[20,] -2.947139 39.326334
[21,] 25.184140 -2.947139
[22,] 11.537432 25.184140
[23,] 213.216987 11.537432
[24,] 11.958219 213.216987
[25,] 159.342378 11.958219
[26,] 38.469194 159.342378
[27,] 30.404768 38.469194
[28,] 5.356646 30.404768
[29,] 153.846869 5.356646
[30,] -826.621428 153.846869
[31,] 11.305477 -826.621428
[32,] -18.722151 11.305477
[33,] -19.334067 -18.722151
[34,] -42.269430 -19.334067
[35,] -97.203245 -42.269430
[36,] -1.719179 -97.203245
[37,] -1.073086 -1.719179
[38,] 25.042512 -1.073086
[39,] 15.078294 25.042512
[40,] 3.074535 15.078294
[41,] -11.273729 3.074535
[42,] 12.377644 -11.273729
[43,] -27.311495 12.377644
[44,] 11.987381 -27.311495
[45,] 31.899100 11.987381
[46,] 172.243632 31.899100
[47,] -35.799408 172.243632
[48,] -18.809349 -35.799408
[49,] 10.454677 -18.809349
[50,] -17.569011 10.454677
[51,] -27.628765 -17.569011
[52,] 197.216666 -27.628765
[53,] 24.691026 197.216666
[54,] 169.249548 24.691026
[55,] 60.618109 169.249548
[56,] -29.116856 60.618109
[57,] -48.521524 -29.116856
[58,] 30.723948 -48.521524
[59,] -12.367897 30.723948
[60,] 26.573202 -12.367897
[61,] -778.294498 26.573202
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 1.779383 15.295541
2 159.456986 1.779383
3 180.628282 159.456986
4 -73.847446 180.628282
5 12.286081 -73.847446
6 -14.800252 12.286081
7 19.741014 -14.800252
8 6.059886 19.741014
9 -26.380852 6.059886
10 51.703124 -26.380852
11 7.374241 51.703124
12 -35.400006 7.374241
13 184.691045 -35.400006
14 -19.399271 184.691045
15 -6.648942 -19.399271
16 18.496160 -6.648942
17 -11.525046 18.496160
18 55.893613 -11.525046
19 39.326334 55.893613
20 -2.947139 39.326334
21 25.184140 -2.947139
22 11.537432 25.184140
23 213.216987 11.537432
24 11.958219 213.216987
25 159.342378 11.958219
26 38.469194 159.342378
27 30.404768 38.469194
28 5.356646 30.404768
29 153.846869 5.356646
30 -826.621428 153.846869
31 11.305477 -826.621428
32 -18.722151 11.305477
33 -19.334067 -18.722151
34 -42.269430 -19.334067
35 -97.203245 -42.269430
36 -1.719179 -97.203245
37 -1.073086 -1.719179
38 25.042512 -1.073086
39 15.078294 25.042512
40 3.074535 15.078294
41 -11.273729 3.074535
42 12.377644 -11.273729
43 -27.311495 12.377644
44 11.987381 -27.311495
45 31.899100 11.987381
46 172.243632 31.899100
47 -35.799408 172.243632
48 -18.809349 -35.799408
49 10.454677 -18.809349
50 -17.569011 10.454677
51 -27.628765 -17.569011
52 197.216666 -27.628765
53 24.691026 197.216666
54 169.249548 24.691026
55 60.618109 169.249548
56 -29.116856 60.618109
57 -48.521524 -29.116856
58 30.723948 -48.521524
59 -12.367897 30.723948
60 26.573202 -12.367897
61 -778.294498 26.573202
> 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/7lnd31292960071.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/8lnd31292960071.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/9lnd31292960071.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/10efu61292960071.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/11hxac1292960071.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/123g9i1292960071.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/13h7781292960071.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/14k85e1292960071.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/155rm21292960071.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/16rr281292960071.tab")
+ }
>
> try(system("convert tmp/1pefu1292960071.ps tmp/1pefu1292960071.png",intern=TRUE))
character(0)
> try(system("convert tmp/2pefu1292960071.ps tmp/2pefu1292960071.png",intern=TRUE))
character(0)
> try(system("convert tmp/3i5ef1292960071.ps tmp/3i5ef1292960071.png",intern=TRUE))
character(0)
> try(system("convert tmp/4i5ef1292960071.ps tmp/4i5ef1292960071.png",intern=TRUE))
character(0)
> try(system("convert tmp/5i5ef1292960071.ps tmp/5i5ef1292960071.png",intern=TRUE))
character(0)
> try(system("convert tmp/6bew01292960071.ps tmp/6bew01292960071.png",intern=TRUE))
character(0)
> try(system("convert tmp/7lnd31292960071.ps tmp/7lnd31292960071.png",intern=TRUE))
character(0)
> try(system("convert tmp/8lnd31292960071.ps tmp/8lnd31292960071.png",intern=TRUE))
character(0)
> try(system("convert tmp/9lnd31292960071.ps tmp/9lnd31292960071.png",intern=TRUE))
character(0)
> try(system("convert tmp/10efu61292960071.ps tmp/10efu61292960071.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
4.029 2.509 4.372