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.0
+ ,-999.0
+ ,38.6
+ ,6654.000
+ ,5712.000
+ ,645.0
+ ,3
+ ,5
+ ,3
+ ,6.3
+ ,2.0
+ ,4.5
+ ,1.000
+ ,6.600
+ ,42.0
+ ,3
+ ,1
+ ,3
+ ,-999.0
+ ,-999.0
+ ,14.0
+ ,3.385
+ ,44.500
+ ,60.0
+ ,1
+ ,1
+ ,1
+ ,-999.0
+ ,-999.0
+ ,-999.0
+ ,0.920
+ ,5.700
+ ,25.0
+ ,5
+ ,2
+ ,3
+ ,2.1
+ ,1.8
+ ,69.0
+ ,2547.000
+ ,4603.000
+ ,624.0
+ ,3
+ ,5
+ ,4
+ ,9.1
+ ,0.7
+ ,27.0
+ ,10.550
+ ,179.500
+ ,180.0
+ ,4
+ ,4
+ ,4
+ ,15.8
+ ,3.9
+ ,19.0
+ ,0.023
+ ,0.300
+ ,35.0
+ ,1
+ ,1
+ ,1
+ ,5.2
+ ,1.0
+ ,30.4
+ ,160.000
+ ,169.000
+ ,392.0
+ ,4
+ ,5
+ ,4
+ ,10.9
+ ,36.0
+ ,28.0
+ ,3.300
+ ,25.600
+ ,63.0
+ ,1
+ ,2
+ ,1
+ ,8.3
+ ,1.4
+ ,50.0
+ ,52.160
+ ,440.000
+ ,230.0
+ ,1
+ ,1
+ ,1
+ ,11.0
+ ,1.5
+ ,7.0
+ ,0.425
+ ,6.400
+ ,112.0
+ ,5
+ ,4
+ ,4
+ ,3.2
+ ,0.7
+ ,30.0
+ ,465.000
+ ,423.000
+ ,281.0
+ ,5
+ ,5
+ ,5
+ ,7.6
+ ,2.7
+ ,-999.0
+ ,0.550
+ ,2.400
+ ,-999.0
+ ,2
+ ,1
+ ,2
+ ,-999.0
+ ,-999.0
+ ,40.0
+ ,187.100
+ ,419.000
+ ,365.0
+ ,5
+ ,5
+ ,5
+ ,6.3
+ ,2.1
+ ,3.5
+ ,0.075
+ ,1.200
+ ,42.0
+ ,1
+ ,1
+ ,1
+ ,8.6
+ ,0.0
+ ,50.0
+ ,3.000
+ ,25.000
+ ,28.0
+ ,2
+ ,2
+ ,2
+ ,6.6
+ ,4.1
+ ,6.0
+ ,0.785
+ ,3.500
+ ,42.0
+ ,2
+ ,2
+ ,2
+ ,9.5
+ ,1.2
+ ,10.4
+ ,0.200
+ ,5.000
+ ,120.0
+ ,2
+ ,2
+ ,2
+ ,4.8
+ ,1.3
+ ,34.0
+ ,1.410
+ ,17.500
+ ,-999.0
+ ,1
+ ,2
+ ,1
+ ,12.0
+ ,6.1
+ ,7.0
+ ,60.000
+ ,81.000
+ ,-999.0
+ ,1
+ ,1
+ ,1
+ ,-999.0
+ ,0.3
+ ,28.0
+ ,529.000
+ ,680.000
+ ,400.0
+ ,5
+ ,5
+ ,5
+ ,3.3
+ ,0.5
+ ,20.0
+ ,27.660
+ ,115.000
+ ,148.0
+ ,5
+ ,5
+ ,5
+ ,11.0
+ ,3.4
+ ,3.9
+ ,0.120
+ ,1.000
+ ,16.0
+ ,3
+ ,1
+ ,2
+ ,-999.0
+ ,-999.0
+ ,39.3
+ ,207.000
+ ,406.000
+ ,252.0
+ ,1
+ ,4
+ ,1
+ ,4.7
+ ,1.5
+ ,41.0
+ ,85.000
+ ,325.000
+ ,310.0
+ ,1
+ ,3
+ ,1
+ ,-999.0
+ ,-999.0
+ ,16.2
+ ,36.330
+ ,119.500
+ ,63.0
+ ,1
+ ,1
+ ,1
+ ,10.4
+ ,3.4
+ ,9.0
+ ,0.101
+ ,4.000
+ ,28.0
+ ,5
+ ,1
+ ,3
+ ,7.4
+ ,0.8
+ ,7.6
+ ,1.040
+ ,5.500
+ ,68.0
+ ,5
+ ,3
+ ,4
+ ,2.1
+ ,0.8
+ ,46.0
+ ,521.000
+ ,655.000
+ ,336.0
+ ,5
+ ,5
+ ,5
+ ,-999.0
+ ,-999.0
+ ,22.4
+ ,100.000
+ ,157.000
+ ,100.0
+ ,1
+ ,1
+ ,1
+ ,-999.0
+ ,-999.0
+ ,16.3
+ ,35.000
+ ,56.000
+ ,33.0
+ ,3
+ ,5
+ ,4
+ ,7.7
+ ,1.4
+ ,2.6
+ ,0.005
+ ,0.140
+ ,21.5
+ ,5
+ ,2
+ ,4
+ ,17.9
+ ,2.0
+ ,24.0
+ ,0.010
+ ,0.250
+ ,50.0
+ ,1
+ ,1
+ ,1
+ ,6.1
+ ,1.9
+ ,100.0
+ ,62.000
+ ,1320.000
+ ,267.0
+ ,1
+ ,1
+ ,1
+ ,8.2
+ ,2.4
+ ,-999.0
+ ,0.122
+ ,3.000
+ ,30.0
+ ,2
+ ,1
+ ,1
+ ,8.4
+ ,2.8
+ ,-999.0
+ ,1.350
+ ,8.100
+ ,45.0
+ ,3
+ ,1
+ ,3
+ ,11.9
+ ,1.3
+ ,3.2
+ ,0.023
+ ,0.400
+ ,19.0
+ ,4
+ ,1
+ ,3
+ ,10.8
+ ,2.0
+ ,2.0
+ ,0.048
+ ,0.330
+ ,30.0
+ ,4
+ ,1
+ ,3
+ ,13.8
+ ,5.6
+ ,5.0
+ ,1.700
+ ,6.300
+ ,12.0
+ ,2
+ ,1
+ ,1
+ ,14.3
+ ,3.1
+ ,6.5
+ ,3.500
+ ,10.800
+ ,120.0
+ ,2
+ ,1
+ ,1
+ ,-999.0
+ ,1.0
+ ,23.6
+ ,250.000
+ ,490.000
+ ,440.0
+ ,5
+ ,5
+ ,5
+ ,15.2
+ ,1.8
+ ,12.0
+ ,0.480
+ ,15.500
+ ,140.0
+ ,2
+ ,2
+ ,2
+ ,10.0
+ ,0.9
+ ,20.2
+ ,10.000
+ ,115.000
+ ,170.0
+ ,4
+ ,4
+ ,4
+ ,11.9
+ ,1.8
+ ,13.0
+ ,1.620
+ ,11.400
+ ,17.0
+ ,2
+ ,1
+ ,2
+ ,6.5
+ ,1.9
+ ,27.0
+ ,192.000
+ ,180.000
+ ,115.0
+ ,4
+ ,4
+ ,4
+ ,7.5
+ ,0.9
+ ,18.0
+ ,2.500
+ ,12.100
+ ,31.0
+ ,5
+ ,5
+ ,5
+ ,-999.0
+ ,-999.0
+ ,13.7
+ ,4.288
+ ,39.200
+ ,63.0
+ ,2
+ ,2
+ ,2
+ ,10.6
+ ,2.6
+ ,4.7
+ ,0.280
+ ,1.900
+ ,21.0
+ ,3
+ ,1
+ ,3
+ ,7.4
+ ,2.4
+ ,9.8
+ ,4.235
+ ,50.400
+ ,52.0
+ ,1
+ ,1
+ ,1
+ ,8.4
+ ,1.2
+ ,29.0
+ ,6.800
+ ,179.000
+ ,164.0
+ ,2
+ ,3
+ ,2
+ ,5.7
+ ,0.9
+ ,7.0
+ ,0.750
+ ,12.300
+ ,225.0
+ ,2
+ ,2
+ ,2
+ ,4.9
+ ,0.5
+ ,6.0
+ ,3.600
+ ,21.000
+ ,225.0
+ ,3
+ ,2
+ ,3
+ ,-999.0
+ ,-999.0
+ ,17.0
+ ,14.830
+ ,98.200
+ ,150.0
+ ,5
+ ,5
+ ,5
+ ,3.2
+ ,0.6
+ ,20.0
+ ,55.500
+ ,175.000
+ ,151.0
+ ,5
+ ,5
+ ,5
+ ,-999.0
+ ,-999.0
+ ,12.7
+ ,1.400
+ ,12.500
+ ,90.0
+ ,2
+ ,2
+ ,2
+ ,8.1
+ ,2.2
+ ,3.5
+ ,0.060
+ ,1.000
+ ,-999.0
+ ,3
+ ,1
+ ,1
+ ,11.0
+ ,2.3
+ ,4.5
+ ,0.900
+ ,2.600
+ ,38.0
+ ,2
+ ,1
+ ,2
+ ,-999.0
+ ,0.5
+ ,7.5
+ ,2.000
+ ,17.000
+ ,200.0
+ ,3
+ ,1
+ ,3
+ ,13.2
+ ,-999.0
+ ,2.3
+ ,0.104
+ ,2.500
+ ,46.0
+ ,3
+ ,2
+ ,2
+ ,9.7
+ ,0.6
+ ,13.0
+ ,4.050
+ ,58.000
+ ,210.0
+ ,4
+ ,3
+ ,4
+ ,12.8
+ ,6.6
+ ,3.0
+ ,3.500
+ ,3.900
+ ,14.0
+ ,2
+ ,1
+ ,1
+ ,4.9
+ ,2.6
+ ,24.0
+ ,4.190
+ ,12.300
+ ,60.0
+ ,3
+ ,1
+ ,2)
+ ,dim=c(9
+ ,62)
+ ,dimnames=list(c('SWS'
+ ,'PS'
+ ,'L'
+ ,'wb'
+ ,'wbr'
+ ,'tg'
+ ,'P'
+ ,'S'
+ ,'D
')
+ ,1:62))
> y <- array(NA,dim=c(9,62),dimnames=list(c('SWS','PS','L','wb','wbr','tg','P','S','D
'),1:62))
> 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 = 'Do not include Seasonal Dummies'
> par1 = '2'
> #'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
PS SWS L wb wbr tg P S D\r
1 -999.0 -999.0 38.6 6654.000 5712.00 645.0 3 5 3
2 2.0 6.3 4.5 1.000 6.60 42.0 3 1 3
3 -999.0 -999.0 14.0 3.385 44.50 60.0 1 1 1
4 -999.0 -999.0 -999.0 0.920 5.70 25.0 5 2 3
5 1.8 2.1 69.0 2547.000 4603.00 624.0 3 5 4
6 0.7 9.1 27.0 10.550 179.50 180.0 4 4 4
7 3.9 15.8 19.0 0.023 0.30 35.0 1 1 1
8 1.0 5.2 30.4 160.000 169.00 392.0 4 5 4
9 36.0 10.9 28.0 3.300 25.60 63.0 1 2 1
10 1.4 8.3 50.0 52.160 440.00 230.0 1 1 1
11 1.5 11.0 7.0 0.425 6.40 112.0 5 4 4
12 0.7 3.2 30.0 465.000 423.00 281.0 5 5 5
13 2.7 7.6 -999.0 0.550 2.40 -999.0 2 1 2
14 -999.0 -999.0 40.0 187.100 419.00 365.0 5 5 5
15 2.1 6.3 3.5 0.075 1.20 42.0 1 1 1
16 0.0 8.6 50.0 3.000 25.00 28.0 2 2 2
17 4.1 6.6 6.0 0.785 3.50 42.0 2 2 2
18 1.2 9.5 10.4 0.200 5.00 120.0 2 2 2
19 1.3 4.8 34.0 1.410 17.50 -999.0 1 2 1
20 6.1 12.0 7.0 60.000 81.00 -999.0 1 1 1
21 0.3 -999.0 28.0 529.000 680.00 400.0 5 5 5
22 0.5 3.3 20.0 27.660 115.00 148.0 5 5 5
23 3.4 11.0 3.9 0.120 1.00 16.0 3 1 2
24 -999.0 -999.0 39.3 207.000 406.00 252.0 1 4 1
25 1.5 4.7 41.0 85.000 325.00 310.0 1 3 1
26 -999.0 -999.0 16.2 36.330 119.50 63.0 1 1 1
27 3.4 10.4 9.0 0.101 4.00 28.0 5 1 3
28 0.8 7.4 7.6 1.040 5.50 68.0 5 3 4
29 0.8 2.1 46.0 521.000 655.00 336.0 5 5 5
30 -999.0 -999.0 22.4 100.000 157.00 100.0 1 1 1
31 -999.0 -999.0 16.3 35.000 56.00 33.0 3 5 4
32 1.4 7.7 2.6 0.005 0.14 21.5 5 2 4
33 2.0 17.9 24.0 0.010 0.25 50.0 1 1 1
34 1.9 6.1 100.0 62.000 1320.00 267.0 1 1 1
35 2.4 8.2 -999.0 0.122 3.00 30.0 2 1 1
36 2.8 8.4 -999.0 1.350 8.10 45.0 3 1 3
37 1.3 11.9 3.2 0.023 0.40 19.0 4 1 3
38 2.0 10.8 2.0 0.048 0.33 30.0 4 1 3
39 5.6 13.8 5.0 1.700 6.30 12.0 2 1 1
40 3.1 14.3 6.5 3.500 10.80 120.0 2 1 1
41 1.0 -999.0 23.6 250.000 490.00 440.0 5 5 5
42 1.8 15.2 12.0 0.480 15.50 140.0 2 2 2
43 0.9 10.0 20.2 10.000 115.00 170.0 4 4 4
44 1.8 11.9 13.0 1.620 11.40 17.0 2 1 2
45 1.9 6.5 27.0 192.000 180.00 115.0 4 4 4
46 0.9 7.5 18.0 2.500 12.10 31.0 5 5 5
47 -999.0 -999.0 13.7 4.288 39.20 63.0 2 2 2
48 2.6 10.6 4.7 0.280 1.90 21.0 3 1 3
49 2.4 7.4 9.8 4.235 50.40 52.0 1 1 1
50 1.2 8.4 29.0 6.800 179.00 164.0 2 3 2
51 0.9 5.7 7.0 0.750 12.30 225.0 2 2 2
52 0.5 4.9 6.0 3.600 21.00 225.0 3 2 3
53 -999.0 -999.0 17.0 14.830 98.20 150.0 5 5 5
54 0.6 3.2 20.0 55.500 175.00 151.0 5 5 5
55 -999.0 -999.0 12.7 1.400 12.50 90.0 2 2 2
56 2.2 8.1 3.5 0.060 1.00 -999.0 3 1 1
57 2.3 11.0 4.5 0.900 2.60 38.0 2 1 2
58 0.5 -999.0 7.5 2.000 17.00 200.0 3 1 3
59 -999.0 13.2 2.3 0.104 2.50 46.0 3 2 2
60 0.6 9.7 13.0 4.050 58.00 210.0 4 3 4
61 6.6 12.8 3.0 3.500 3.90 14.0 2 1 1
62 2.6 4.9 24.0 4.190 12.30 60.0 3 1 2
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) SWS L wb wbr tg
-90.747120 0.747139 0.048076 -0.058508 0.050150 0.007202
P S `D\r`
-56.503995 -57.168937 138.812334
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-912.45 -53.93 -14.48 55.57 704.36
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -90.747120 72.068777 -1.259 0.2135
SWS 0.747139 0.077462 9.645 2.93e-13 ***
L 0.048076 0.129588 0.371 0.7121
wb -0.058508 0.099518 -0.588 0.5591
wbr 0.050150 0.099368 0.505 0.6159
tg 0.007202 0.118331 0.061 0.9517
P -56.503995 59.041082 -0.957 0.3429
S -57.168937 37.237168 -1.535 0.1307
`D\r` 138.812334 75.848946 1.830 0.0729 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 236.7 on 53 degrees of freedom
Multiple R-squared: 0.6943, Adjusted R-squared: 0.6481
F-statistic: 15.05 on 8 and 53 DF, p-value: 3.107e-11
> 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,] 2.469785e-04 4.939569e-04 0.9997530
[2,] 1.993374e-05 3.986749e-05 0.9999801
[3,] 1.088382e-06 2.176765e-06 0.9999989
[4,] 5.024930e-08 1.004986e-07 0.9999999
[5,] 2.355591e-09 4.711183e-09 1.0000000
[6,] 9.158024e-11 1.831605e-10 1.0000000
[7,] 3.950450e-12 7.900900e-12 1.0000000
[8,] 1.476013e-13 2.952025e-13 1.0000000
[9,] 5.061939e-15 1.012388e-14 1.0000000
[10,] 3.630902e-01 7.261803e-01 0.6369098
[11,] 2.900966e-01 5.801932e-01 0.7099034
[12,] 2.215534e-01 4.431067e-01 0.7784466
[13,] 1.676603e-01 3.353205e-01 0.8323397
[14,] 1.340730e-01 2.681459e-01 0.8659270
[15,] 1.126999e-01 2.253997e-01 0.8873001
[16,] 7.632347e-02 1.526469e-01 0.9236765
[17,] 4.945799e-02 9.891599e-02 0.9505420
[18,] 5.150218e-02 1.030044e-01 0.9484978
[19,] 5.909410e-02 1.181882e-01 0.9409059
[20,] 5.741794e-02 1.148359e-01 0.9425821
[21,] 3.706056e-02 7.412112e-02 0.9629394
[22,] 2.301069e-02 4.602138e-02 0.9769893
[23,] 3.432208e-02 6.864417e-02 0.9656779
[24,] 3.102232e-02 6.204464e-02 0.9689777
[25,] 1.984750e-02 3.969501e-02 0.9801525
[26,] 1.207691e-02 2.415382e-02 0.9879231
[27,] 7.397858e-03 1.479572e-02 0.9926021
[28,] 5.218207e-03 1.043641e-02 0.9947818
[29,] 4.675678e-03 9.351355e-03 0.9953243
[30,] 5.710477e-02 1.142095e-01 0.9428952
[31,] 3.920766e-02 7.841531e-02 0.9607923
[32,] 2.357352e-02 4.714704e-02 0.9764265
[33,] 1.590174e-02 3.180348e-02 0.9840983
[34,] 8.767056e-03 1.753411e-02 0.9912329
[35,] 1.831791e-02 3.663583e-02 0.9816821
[36,] 1.799898e-02 3.599796e-02 0.9820010
[37,] 2.623105e-02 5.246210e-02 0.9737689
[38,] 1.553613e-02 3.107225e-02 0.9844639
[39,] 5.841097e-02 1.168219e-01 0.9415890
> postscript(file="/var/www/html/freestat/rcomp/tmp/1beia1292871063.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/2beia1292871063.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/3beia1292871063.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/445zd1292871063.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/545zd1292871063.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
-26.582934 -102.507232 -190.138963 -133.824473 -98.544272 -26.888435
7 8 9 10 11 12
56.523720 41.074754 147.742252 37.732174 38.535755 -34.113675
13 14 15 16 17 18
35.455051 -302.175305 62.474210 30.359270 38.916674 32.967247
19 20 21 22 23 24
125.255270 69.048466 704.364453 -43.091706 34.638991 -27.447170
25 26 27 28 29 30
162.404923 -192.100010 8.799763 -16.274330 -42.715407 -190.819934
31 32 33 34 35 36
-263.535334 -72.283902 52.708070 -6.350431 166.053339 -55.108704
37 38 39 40 41 42
-50.405311 -48.900014 116.857913 113.014079 698.192543 28.577406
43 44 45 46 47 48
-23.759449 -25.015915 -12.686495 -41.202591 -214.966921 -104.784731
49 50 51 52 53 54
59.353489 81.407049 34.579732 -47.752374 -293.512372 -44.318711
55 56 57 58 59 60
-213.943267 181.743555 -23.186892 645.346686 -912.451072 -78.435759
61 62
118.912475 36.784779
> postscript(file="/var/www/html/freestat/rcomp/tmp/645zd1292871063.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 -26.582934 NA
1 -102.507232 -26.582934
2 -190.138963 -102.507232
3 -133.824473 -190.138963
4 -98.544272 -133.824473
5 -26.888435 -98.544272
6 56.523720 -26.888435
7 41.074754 56.523720
8 147.742252 41.074754
9 37.732174 147.742252
10 38.535755 37.732174
11 -34.113675 38.535755
12 35.455051 -34.113675
13 -302.175305 35.455051
14 62.474210 -302.175305
15 30.359270 62.474210
16 38.916674 30.359270
17 32.967247 38.916674
18 125.255270 32.967247
19 69.048466 125.255270
20 704.364453 69.048466
21 -43.091706 704.364453
22 34.638991 -43.091706
23 -27.447170 34.638991
24 162.404923 -27.447170
25 -192.100010 162.404923
26 8.799763 -192.100010
27 -16.274330 8.799763
28 -42.715407 -16.274330
29 -190.819934 -42.715407
30 -263.535334 -190.819934
31 -72.283902 -263.535334
32 52.708070 -72.283902
33 -6.350431 52.708070
34 166.053339 -6.350431
35 -55.108704 166.053339
36 -50.405311 -55.108704
37 -48.900014 -50.405311
38 116.857913 -48.900014
39 113.014079 116.857913
40 698.192543 113.014079
41 28.577406 698.192543
42 -23.759449 28.577406
43 -25.015915 -23.759449
44 -12.686495 -25.015915
45 -41.202591 -12.686495
46 -214.966921 -41.202591
47 -104.784731 -214.966921
48 59.353489 -104.784731
49 81.407049 59.353489
50 34.579732 81.407049
51 -47.752374 34.579732
52 -293.512372 -47.752374
53 -44.318711 -293.512372
54 -213.943267 -44.318711
55 181.743555 -213.943267
56 -23.186892 181.743555
57 645.346686 -23.186892
58 -912.451072 645.346686
59 -78.435759 -912.451072
60 118.912475 -78.435759
61 36.784779 118.912475
62 NA 36.784779
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -102.507232 -26.582934
[2,] -190.138963 -102.507232
[3,] -133.824473 -190.138963
[4,] -98.544272 -133.824473
[5,] -26.888435 -98.544272
[6,] 56.523720 -26.888435
[7,] 41.074754 56.523720
[8,] 147.742252 41.074754
[9,] 37.732174 147.742252
[10,] 38.535755 37.732174
[11,] -34.113675 38.535755
[12,] 35.455051 -34.113675
[13,] -302.175305 35.455051
[14,] 62.474210 -302.175305
[15,] 30.359270 62.474210
[16,] 38.916674 30.359270
[17,] 32.967247 38.916674
[18,] 125.255270 32.967247
[19,] 69.048466 125.255270
[20,] 704.364453 69.048466
[21,] -43.091706 704.364453
[22,] 34.638991 -43.091706
[23,] -27.447170 34.638991
[24,] 162.404923 -27.447170
[25,] -192.100010 162.404923
[26,] 8.799763 -192.100010
[27,] -16.274330 8.799763
[28,] -42.715407 -16.274330
[29,] -190.819934 -42.715407
[30,] -263.535334 -190.819934
[31,] -72.283902 -263.535334
[32,] 52.708070 -72.283902
[33,] -6.350431 52.708070
[34,] 166.053339 -6.350431
[35,] -55.108704 166.053339
[36,] -50.405311 -55.108704
[37,] -48.900014 -50.405311
[38,] 116.857913 -48.900014
[39,] 113.014079 116.857913
[40,] 698.192543 113.014079
[41,] 28.577406 698.192543
[42,] -23.759449 28.577406
[43,] -25.015915 -23.759449
[44,] -12.686495 -25.015915
[45,] -41.202591 -12.686495
[46,] -214.966921 -41.202591
[47,] -104.784731 -214.966921
[48,] 59.353489 -104.784731
[49,] 81.407049 59.353489
[50,] 34.579732 81.407049
[51,] -47.752374 34.579732
[52,] -293.512372 -47.752374
[53,] -44.318711 -293.512372
[54,] -213.943267 -44.318711
[55,] 181.743555 -213.943267
[56,] -23.186892 181.743555
[57,] 645.346686 -23.186892
[58,] -912.451072 645.346686
[59,] -78.435759 -912.451072
[60,] 118.912475 -78.435759
[61,] 36.784779 118.912475
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -102.507232 -26.582934
2 -190.138963 -102.507232
3 -133.824473 -190.138963
4 -98.544272 -133.824473
5 -26.888435 -98.544272
6 56.523720 -26.888435
7 41.074754 56.523720
8 147.742252 41.074754
9 37.732174 147.742252
10 38.535755 37.732174
11 -34.113675 38.535755
12 35.455051 -34.113675
13 -302.175305 35.455051
14 62.474210 -302.175305
15 30.359270 62.474210
16 38.916674 30.359270
17 32.967247 38.916674
18 125.255270 32.967247
19 69.048466 125.255270
20 704.364453 69.048466
21 -43.091706 704.364453
22 34.638991 -43.091706
23 -27.447170 34.638991
24 162.404923 -27.447170
25 -192.100010 162.404923
26 8.799763 -192.100010
27 -16.274330 8.799763
28 -42.715407 -16.274330
29 -190.819934 -42.715407
30 -263.535334 -190.819934
31 -72.283902 -263.535334
32 52.708070 -72.283902
33 -6.350431 52.708070
34 166.053339 -6.350431
35 -55.108704 166.053339
36 -50.405311 -55.108704
37 -48.900014 -50.405311
38 116.857913 -48.900014
39 113.014079 116.857913
40 698.192543 113.014079
41 28.577406 698.192543
42 -23.759449 28.577406
43 -25.015915 -23.759449
44 -12.686495 -25.015915
45 -41.202591 -12.686495
46 -214.966921 -41.202591
47 -104.784731 -214.966921
48 59.353489 -104.784731
49 81.407049 59.353489
50 34.579732 81.407049
51 -47.752374 34.579732
52 -293.512372 -47.752374
53 -44.318711 -293.512372
54 -213.943267 -44.318711
55 181.743555 -213.943267
56 -23.186892 181.743555
57 645.346686 -23.186892
58 -912.451072 645.346686
59 -78.435759 -912.451072
60 118.912475 -78.435759
61 36.784779 118.912475
> 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/7wfyy1292871063.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/8wfyy1292871063.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/976x11292871063.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')
Warning messages:
1: In sqrt(crit * p * (1 - hh)/hh) : NaNs produced
2: In sqrt(crit * p * (1 - hh)/hh) : NaNs produced
> par(opar)
> dev.off()
null device
1
> if (n > n25) {
+ postscript(file="/var/www/html/freestat/rcomp/tmp/1076x11292871063.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/11s6e71292871063.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/12epcd1292871063.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/13lq961292871063.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/14vh9r1292871063.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/15zipx1292871063.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/16vs561292871063.tab")
+ }
>
> try(system("convert tmp/1beia1292871063.ps tmp/1beia1292871063.png",intern=TRUE))
character(0)
> try(system("convert tmp/2beia1292871063.ps tmp/2beia1292871063.png",intern=TRUE))
character(0)
> try(system("convert tmp/3beia1292871063.ps tmp/3beia1292871063.png",intern=TRUE))
character(0)
> try(system("convert tmp/445zd1292871063.ps tmp/445zd1292871063.png",intern=TRUE))
character(0)
> try(system("convert tmp/545zd1292871063.ps tmp/545zd1292871063.png",intern=TRUE))
character(0)
> try(system("convert tmp/645zd1292871063.ps tmp/645zd1292871063.png",intern=TRUE))
character(0)
> try(system("convert tmp/7wfyy1292871063.ps tmp/7wfyy1292871063.png",intern=TRUE))
character(0)
> try(system("convert tmp/8wfyy1292871063.ps tmp/8wfyy1292871063.png",intern=TRUE))
character(0)
> try(system("convert tmp/976x11292871063.ps tmp/976x11292871063.png",intern=TRUE))
character(0)
> try(system("convert tmp/1076x11292871063.ps tmp/1076x11292871063.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.895 2.420 4.490