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(19435.1
+ ,2.01
+ ,20604.6
+ ,20604.6
+ ,20604.6
+ ,22686.8
+ ,2.01
+ ,19435.1
+ ,18714.9
+ ,18714.9
+ ,20396.7
+ ,2.01
+ ,22686.8
+ ,19435.1
+ ,18492.6
+ ,19233.6
+ ,2.01
+ ,20396.7
+ ,22686.8
+ ,19435.1
+ ,22751
+ ,2.01
+ ,19233.6
+ ,20396.7
+ ,22686.8
+ ,19864
+ ,2.01
+ ,22751
+ ,19233.6
+ ,20396.7
+ ,17165.4
+ ,2.02
+ ,19864
+ ,22751
+ ,19233.6
+ ,22309.7
+ ,2.02
+ ,17165.4
+ ,19864
+ ,22751
+ ,21786.3
+ ,2.03
+ ,22309.7
+ ,17165.4
+ ,19864
+ ,21927.6
+ ,2.05
+ ,21786.3
+ ,22309.7
+ ,17165.4
+ ,20957.9
+ ,2.08
+ ,21927.6
+ ,21786.3
+ ,22309.7
+ ,19726
+ ,2.07
+ ,20957.9
+ ,21927.6
+ ,21786.3
+ ,21315.7
+ ,2.06
+ ,19726
+ ,20957.9
+ ,21927.6
+ ,24771.5
+ ,2.05
+ ,21315.7
+ ,19726
+ ,20957.9
+ ,22592.4
+ ,2.05
+ ,24771.5
+ ,21315.7
+ ,19726
+ ,21942.1
+ ,2.05
+ ,22592.4
+ ,24771.5
+ ,21315.7
+ ,23973.7
+ ,2.05
+ ,21942.1
+ ,22592.4
+ ,24771.5
+ ,20815.7
+ ,2.05
+ ,23973.7
+ ,21942.1
+ ,22592.4
+ ,19931.4
+ ,2.06
+ ,20815.7
+ ,23973.7
+ ,21942.1
+ ,24436.8
+ ,2.06
+ ,19931.4
+ ,20815.7
+ ,23973.7
+ ,22838.7
+ ,2.07
+ ,24436.8
+ ,19931.4
+ ,20815.7
+ ,24465.3
+ ,2.07
+ ,22838.7
+ ,24436.8
+ ,19931.4
+ ,23007.3
+ ,2.3
+ ,24465.3
+ ,22838.7
+ ,24436.8
+ ,22720.8
+ ,2.31
+ ,23007.3
+ ,24465.3
+ ,22838.7
+ ,23045.7
+ ,2.31
+ ,22720.8
+ ,23007.3
+ ,24465.3
+ ,27198.5
+ ,2.53
+ ,23045.7
+ ,22720.8
+ ,23007.3
+ ,22401.9
+ ,2.58
+ ,27198.5
+ ,23045.7
+ ,22720.8
+ ,25122.7
+ ,2.59
+ ,22401.9
+ ,27198.5
+ ,23045.7
+ ,26100.5
+ ,2.73
+ ,25122.7
+ ,22401.9
+ ,27198.5
+ ,22904.9
+ ,2.82
+ ,26100.5
+ ,25122.7
+ ,22401.9
+ ,22040.4
+ ,3
+ ,22904.9
+ ,26100.5
+ ,25122.7
+ ,25981.5
+ ,3.04
+ ,22040.4
+ ,22904.9
+ ,26100.5
+ ,26157.1
+ ,3.23
+ ,25981.5
+ ,22040.4
+ ,22904.9
+ ,25975.4
+ ,3.32
+ ,26157.1
+ ,25981.5
+ ,22040.4
+ ,22589.8
+ ,3.49
+ ,25975.4
+ ,26157.1
+ ,25981.5
+ ,25370.4
+ ,3.57
+ ,22589.8
+ ,25975.4
+ ,26157.1
+ ,25091.1
+ ,3.56
+ ,25370.4
+ ,22589.8
+ ,25975.4
+ ,28760.9
+ ,3.72
+ ,25091.1
+ ,25370.4
+ ,22589.8
+ ,24325.9
+ ,3.82
+ ,28760.9
+ ,25091.1
+ ,25370.4
+ ,25821.7
+ ,3.82
+ ,24325.9
+ ,28760.9
+ ,25091.1
+ ,27645.7
+ ,3.98
+ ,25821.7
+ ,24325.9
+ ,28760.9
+ ,26296.9
+ ,4.06
+ ,27645.7
+ ,25821.7
+ ,24325.9
+ ,24141.5
+ ,4.08
+ ,26296.9
+ ,27645.7
+ ,25821.7
+ ,27268.1
+ ,4.19
+ ,24141.5
+ ,26296.9
+ ,27645.7
+ ,29060.3
+ ,4.16
+ ,27268.1
+ ,24141.5
+ ,26296.9
+ ,28226.4
+ ,4.17
+ ,29060.3
+ ,27268.1
+ ,24141.5
+ ,23268.5
+ ,4.21
+ ,28226.4
+ ,29060.3
+ ,27268.1
+ ,26938.2
+ ,4.21
+ ,23268.5
+ ,28226.4
+ ,29060.3
+ ,27217.5
+ ,4.17
+ ,26938.2
+ ,23268.5
+ ,28226.4
+ ,27540.5
+ ,4.19
+ ,27217.5
+ ,26938.2
+ ,23268.5
+ ,29167.6
+ ,4.25
+ ,27540.5
+ ,27217.5
+ ,26938.2
+ ,26671.5
+ ,4.25
+ ,29167.6
+ ,27540.5
+ ,27217.5
+ ,30184
+ ,4.2
+ ,26671.5
+ ,29167.6
+ ,27540.5
+ ,28422.3
+ ,4.33
+ ,30184
+ ,26671.5
+ ,29167.6
+ ,23774.3
+ ,4.41
+ ,28422.3
+ ,30184
+ ,26671.5
+ ,29601
+ ,4.56
+ ,23774.3
+ ,28422.3
+ ,30184
+ ,28523.6
+ ,5.18
+ ,29601
+ ,23774.3
+ ,28422.3
+ ,23622
+ ,3.42
+ ,28523.6
+ ,29601
+ ,23774.3
+ ,21320.3
+ ,2.71
+ ,23622
+ ,28523.6
+ ,29601
+ ,20423.6
+ ,2.29
+ ,21320.3
+ ,23622
+ ,28523.6
+ ,21174.9
+ ,2
+ ,20423.6
+ ,21320.3
+ ,23622
+ ,23050.2
+ ,1.64
+ ,21174.9
+ ,20423.6
+ ,21320.3
+ ,21202.9
+ ,1.3
+ ,23050.2
+ ,21174.9
+ ,20423.6
+ ,20476.4
+ ,1.08
+ ,21202.9
+ ,23050.2
+ ,21174.9
+ ,23173.3
+ ,1
+ ,20476.4
+ ,21202.9
+ ,23050.2
+ ,22468
+ ,1
+ ,23173.3
+ ,20476.4
+ ,21202.9
+ ,19842.7
+ ,1
+ ,22468
+ ,23173.3
+ ,20476.4)
+ ,dim=c(5
+ ,67)
+ ,dimnames=list(c('Y'
+ ,'X'
+ ,'Y1'
+ ,'Y2'
+ ,'Y3')
+ ,1:67))
> y <- array(NA,dim=c(5,67),dimnames=list(c('Y','X','Y1','Y2','Y3'),1:67))
> 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 Y1 Y2 Y3 M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 19435.1 2.01 20604.6 20604.6 20604.6 1 0 0 0 0 0 0 0 0 0 0 1
2 22686.8 2.01 19435.1 18714.9 18714.9 0 1 0 0 0 0 0 0 0 0 0 2
3 20396.7 2.01 22686.8 19435.1 18492.6 0 0 1 0 0 0 0 0 0 0 0 3
4 19233.6 2.01 20396.7 22686.8 19435.1 0 0 0 1 0 0 0 0 0 0 0 4
5 22751.0 2.01 19233.6 20396.7 22686.8 0 0 0 0 1 0 0 0 0 0 0 5
6 19864.0 2.01 22751.0 19233.6 20396.7 0 0 0 0 0 1 0 0 0 0 0 6
7 17165.4 2.02 19864.0 22751.0 19233.6 0 0 0 0 0 0 1 0 0 0 0 7
8 22309.7 2.02 17165.4 19864.0 22751.0 0 0 0 0 0 0 0 1 0 0 0 8
9 21786.3 2.03 22309.7 17165.4 19864.0 0 0 0 0 0 0 0 0 1 0 0 9
10 21927.6 2.05 21786.3 22309.7 17165.4 0 0 0 0 0 0 0 0 0 1 0 10
11 20957.9 2.08 21927.6 21786.3 22309.7 0 0 0 0 0 0 0 0 0 0 1 11
12 19726.0 2.07 20957.9 21927.6 21786.3 0 0 0 0 0 0 0 0 0 0 0 12
13 21315.7 2.06 19726.0 20957.9 21927.6 1 0 0 0 0 0 0 0 0 0 0 13
14 24771.5 2.05 21315.7 19726.0 20957.9 0 1 0 0 0 0 0 0 0 0 0 14
15 22592.4 2.05 24771.5 21315.7 19726.0 0 0 1 0 0 0 0 0 0 0 0 15
16 21942.1 2.05 22592.4 24771.5 21315.7 0 0 0 1 0 0 0 0 0 0 0 16
17 23973.7 2.05 21942.1 22592.4 24771.5 0 0 0 0 1 0 0 0 0 0 0 17
18 20815.7 2.05 23973.7 21942.1 22592.4 0 0 0 0 0 1 0 0 0 0 0 18
19 19931.4 2.06 20815.7 23973.7 21942.1 0 0 0 0 0 0 1 0 0 0 0 19
20 24436.8 2.06 19931.4 20815.7 23973.7 0 0 0 0 0 0 0 1 0 0 0 20
21 22838.7 2.07 24436.8 19931.4 20815.7 0 0 0 0 0 0 0 0 1 0 0 21
22 24465.3 2.07 22838.7 24436.8 19931.4 0 0 0 0 0 0 0 0 0 1 0 22
23 23007.3 2.30 24465.3 22838.7 24436.8 0 0 0 0 0 0 0 0 0 0 1 23
24 22720.8 2.31 23007.3 24465.3 22838.7 0 0 0 0 0 0 0 0 0 0 0 24
25 23045.7 2.31 22720.8 23007.3 24465.3 1 0 0 0 0 0 0 0 0 0 0 25
26 27198.5 2.53 23045.7 22720.8 23007.3 0 1 0 0 0 0 0 0 0 0 0 26
27 22401.9 2.58 27198.5 23045.7 22720.8 0 0 1 0 0 0 0 0 0 0 0 27
28 25122.7 2.59 22401.9 27198.5 23045.7 0 0 0 1 0 0 0 0 0 0 0 28
29 26100.5 2.73 25122.7 22401.9 27198.5 0 0 0 0 1 0 0 0 0 0 0 29
30 22904.9 2.82 26100.5 25122.7 22401.9 0 0 0 0 0 1 0 0 0 0 0 30
31 22040.4 3.00 22904.9 26100.5 25122.7 0 0 0 0 0 0 1 0 0 0 0 31
32 25981.5 3.04 22040.4 22904.9 26100.5 0 0 0 0 0 0 0 1 0 0 0 32
33 26157.1 3.23 25981.5 22040.4 22904.9 0 0 0 0 0 0 0 0 1 0 0 33
34 25975.4 3.32 26157.1 25981.5 22040.4 0 0 0 0 0 0 0 0 0 1 0 34
35 22589.8 3.49 25975.4 26157.1 25981.5 0 0 0 0 0 0 0 0 0 0 1 35
36 25370.4 3.57 22589.8 25975.4 26157.1 0 0 0 0 0 0 0 0 0 0 0 36
37 25091.1 3.56 25370.4 22589.8 25975.4 1 0 0 0 0 0 0 0 0 0 0 37
38 28760.9 3.72 25091.1 25370.4 22589.8 0 1 0 0 0 0 0 0 0 0 0 38
39 24325.9 3.82 28760.9 25091.1 25370.4 0 0 1 0 0 0 0 0 0 0 0 39
40 25821.7 3.82 24325.9 28760.9 25091.1 0 0 0 1 0 0 0 0 0 0 0 40
41 27645.7 3.98 25821.7 24325.9 28760.9 0 0 0 0 1 0 0 0 0 0 0 41
42 26296.9 4.06 27645.7 25821.7 24325.9 0 0 0 0 0 1 0 0 0 0 0 42
43 24141.5 4.08 26296.9 27645.7 25821.7 0 0 0 0 0 0 1 0 0 0 0 43
44 27268.1 4.19 24141.5 26296.9 27645.7 0 0 0 0 0 0 0 1 0 0 0 44
45 29060.3 4.16 27268.1 24141.5 26296.9 0 0 0 0 0 0 0 0 1 0 0 45
46 28226.4 4.17 29060.3 27268.1 24141.5 0 0 0 0 0 0 0 0 0 1 0 46
47 23268.5 4.21 28226.4 29060.3 27268.1 0 0 0 0 0 0 0 0 0 0 1 47
48 26938.2 4.21 23268.5 28226.4 29060.3 0 0 0 0 0 0 0 0 0 0 0 48
49 27217.5 4.17 26938.2 23268.5 28226.4 1 0 0 0 0 0 0 0 0 0 0 49
50 27540.5 4.19 27217.5 26938.2 23268.5 0 1 0 0 0 0 0 0 0 0 0 50
51 29167.6 4.25 27540.5 27217.5 26938.2 0 0 1 0 0 0 0 0 0 0 0 51
52 26671.5 4.25 29167.6 27540.5 27217.5 0 0 0 1 0 0 0 0 0 0 0 52
53 30184.0 4.20 26671.5 29167.6 27540.5 0 0 0 0 1 0 0 0 0 0 0 53
54 28422.3 4.33 30184.0 26671.5 29167.6 0 0 0 0 0 1 0 0 0 0 0 54
55 23774.3 4.41 28422.3 30184.0 26671.5 0 0 0 0 0 0 1 0 0 0 0 55
56 29601.0 4.56 23774.3 28422.3 30184.0 0 0 0 0 0 0 0 1 0 0 0 56
57 28523.6 5.18 29601.0 23774.3 28422.3 0 0 0 0 0 0 0 0 1 0 0 57
58 23622.0 3.42 28523.6 29601.0 23774.3 0 0 0 0 0 0 0 0 0 1 0 58
59 21320.3 2.71 23622.0 28523.6 29601.0 0 0 0 0 0 0 0 0 0 0 1 59
60 20423.6 2.29 21320.3 23622.0 28523.6 0 0 0 0 0 0 0 0 0 0 0 60
61 21174.9 2.00 20423.6 21320.3 23622.0 1 0 0 0 0 0 0 0 0 0 0 61
62 23050.2 1.64 21174.9 20423.6 21320.3 0 1 0 0 0 0 0 0 0 0 0 62
63 21202.9 1.30 23050.2 21174.9 20423.6 0 0 1 0 0 0 0 0 0 0 0 63
64 20476.4 1.08 21202.9 23050.2 21174.9 0 0 0 1 0 0 0 0 0 0 0 64
65 23173.3 1.00 20476.4 21202.9 23050.2 0 0 0 0 1 0 0 0 0 0 0 65
66 22468.0 1.00 23173.3 20476.4 21202.9 0 0 0 0 0 1 0 0 0 0 0 66
67 19842.7 1.00 22468.0 23173.3 20476.4 0 0 0 0 0 0 1 0 0 0 0 67
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X Y1 Y2 Y3 M1
9374.5747 1224.0078 0.0977 0.1160 0.1825 724.3543
M2 M3 M4 M5 M6 M7
3884.9272 1126.7366 815.9470 2970.7550 924.9910 -1503.4400
M8 M9 M10 M11 t
2712.3348 2530.8053 2003.6607 -1276.1919 10.7061
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-3121.8 -838.3 100.8 878.4 2348.6
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 9.375e+03 4.383e+03 2.139 0.037341 *
X 1.224e+03 5.081e+02 2.409 0.019729 *
Y1 9.771e-02 1.551e-01 0.630 0.531598
Y2 1.159e-01 1.483e-01 0.782 0.438003
Y3 1.825e-01 1.525e-01 1.197 0.236961
M1 7.244e+02 8.855e+02 0.818 0.417241
M2 3.885e+03 1.010e+03 3.848 0.000339 ***
M3 1.127e+03 1.124e+03 1.003 0.320766
M4 8.159e+02 9.295e+02 0.878 0.384214
M5 2.971e+03 8.411e+02 3.532 0.000897 ***
M6 9.250e+02 1.056e+03 0.876 0.385160
M7 -1.503e+03 9.105e+02 -1.651 0.104959
M8 2.712e+03 8.806e+02 3.080 0.003359 **
M9 2.531e+03 1.204e+03 2.101 0.040665 *
M10 2.004e+03 1.199e+03 1.671 0.101069
M11 -1.276e+03 9.169e+02 -1.392 0.170113
t 1.071e+01 1.389e+01 0.771 0.444501
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1304 on 50 degrees of freedom
Multiple R-squared: 0.8566, Adjusted R-squared: 0.8107
F-statistic: 18.67 on 16 and 50 DF, p-value: 9.752e-16
> 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.0816056499 0.163211300 0.9183944
[2,] 0.0299617613 0.059923523 0.9700382
[3,] 0.0113498581 0.022699716 0.9886501
[4,] 0.0072290631 0.014458126 0.9927709
[5,] 0.0112282167 0.022456433 0.9887718
[6,] 0.0076527463 0.015305493 0.9923473
[7,] 0.0036915806 0.007383161 0.9963084
[8,] 0.0202067771 0.040413554 0.9797932
[9,] 0.0213375190 0.042675038 0.9786625
[10,] 0.0119109742 0.023821948 0.9880890
[11,] 0.0112884877 0.022576975 0.9887115
[12,] 0.0062173024 0.012434605 0.9937827
[13,] 0.0028930620 0.005786124 0.9971069
[14,] 0.0015115019 0.003023004 0.9984885
[15,] 0.0006533759 0.001306752 0.9993466
[16,] 0.0035927075 0.007185415 0.9964073
[17,] 0.0020926991 0.004185398 0.9979073
[18,] 0.0009180905 0.001836181 0.9990819
[19,] 0.0007463948 0.001492790 0.9992536
[20,] 0.0012892926 0.002578585 0.9987107
[21,] 0.0005826279 0.001165256 0.9994174
[22,] 0.0005595974 0.001119195 0.9994404
[23,] 0.0014995264 0.002999053 0.9985005
[24,] 0.0054343540 0.010868708 0.9945656
[25,] 0.0241537737 0.048307547 0.9758462
[26,] 0.1363072004 0.272614401 0.8636928
[27,] 0.1023433993 0.204686799 0.8976566
[28,] 0.1183131663 0.236626333 0.8816868
> postscript(file="/var/www/html/rcomp/tmp/1ty3q1258477215.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/2sf9t1258477215.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/39qno1258477215.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/4c5mh1258477215.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/5n9r41258477215.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 = 67
Frequency = 1
1 2 3 4 5
-1297.9066872 -539.1936842 -442.4554092 -1630.7839588 -493.2138532
6 7 8 9 10
-1135.9733763 -1342.5736137 -468.3258422 -495.9263473 84.5248815
11 12 13 14 15
1455.1956751 -877.4694952 196.4237011 657.6933635 928.9413763
16 17 18 19 20
100.7743347 -347.6824620 -1195.9904081 516.8739375 877.5651792
21 22 23 24 25
-323.2209084 1614.9509795 2348.6238067 1008.5183688 498.5223432
26 27 28 29 30
1478.3519601 -1023.0993050 1913.3722253 86.6665761 -719.5807012
31 32 33 34 35
315.5728749 257.7689434 670.0639648 578.2918592 -468.1818801
36 37 38 39 40
1247.4117546 399.3430621 1024.8346851 -1618.7742441 235.8907471
41 42 43 44 45
-603.1761451 442.9804375 328.0995669 -872.3461908 1318.0171882
46 47 48 49 50
844.0728174 -1590.6423445 1046.1574729 1007.8906812 -1412.7522133
51 52 53 54 55
2154.6523418 -288.7727098 1115.6772577 879.1719557 -1228.5899052
56 57 58 59 60
205.3379104 -1168.9338974 -3121.8405376 -1744.9952573 -2424.6181011
61 62 63 64 65
-804.2731005 -1208.9341113 0.7352401 -330.4806386 241.7286265
66 67
1729.3920925 1410.6171396
> postscript(file="/var/www/html/rcomp/tmp/6v3f11258477215.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 = 67
Frequency = 1
lag(myerror, k = 1) myerror
0 -1297.9066872 NA
1 -539.1936842 -1297.9066872
2 -442.4554092 -539.1936842
3 -1630.7839588 -442.4554092
4 -493.2138532 -1630.7839588
5 -1135.9733763 -493.2138532
6 -1342.5736137 -1135.9733763
7 -468.3258422 -1342.5736137
8 -495.9263473 -468.3258422
9 84.5248815 -495.9263473
10 1455.1956751 84.5248815
11 -877.4694952 1455.1956751
12 196.4237011 -877.4694952
13 657.6933635 196.4237011
14 928.9413763 657.6933635
15 100.7743347 928.9413763
16 -347.6824620 100.7743347
17 -1195.9904081 -347.6824620
18 516.8739375 -1195.9904081
19 877.5651792 516.8739375
20 -323.2209084 877.5651792
21 1614.9509795 -323.2209084
22 2348.6238067 1614.9509795
23 1008.5183688 2348.6238067
24 498.5223432 1008.5183688
25 1478.3519601 498.5223432
26 -1023.0993050 1478.3519601
27 1913.3722253 -1023.0993050
28 86.6665761 1913.3722253
29 -719.5807012 86.6665761
30 315.5728749 -719.5807012
31 257.7689434 315.5728749
32 670.0639648 257.7689434
33 578.2918592 670.0639648
34 -468.1818801 578.2918592
35 1247.4117546 -468.1818801
36 399.3430621 1247.4117546
37 1024.8346851 399.3430621
38 -1618.7742441 1024.8346851
39 235.8907471 -1618.7742441
40 -603.1761451 235.8907471
41 442.9804375 -603.1761451
42 328.0995669 442.9804375
43 -872.3461908 328.0995669
44 1318.0171882 -872.3461908
45 844.0728174 1318.0171882
46 -1590.6423445 844.0728174
47 1046.1574729 -1590.6423445
48 1007.8906812 1046.1574729
49 -1412.7522133 1007.8906812
50 2154.6523418 -1412.7522133
51 -288.7727098 2154.6523418
52 1115.6772577 -288.7727098
53 879.1719557 1115.6772577
54 -1228.5899052 879.1719557
55 205.3379104 -1228.5899052
56 -1168.9338974 205.3379104
57 -3121.8405376 -1168.9338974
58 -1744.9952573 -3121.8405376
59 -2424.6181011 -1744.9952573
60 -804.2731005 -2424.6181011
61 -1208.9341113 -804.2731005
62 0.7352401 -1208.9341113
63 -330.4806386 0.7352401
64 241.7286265 -330.4806386
65 1729.3920925 241.7286265
66 1410.6171396 1729.3920925
67 NA 1410.6171396
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -539.1936842 -1297.9066872
[2,] -442.4554092 -539.1936842
[3,] -1630.7839588 -442.4554092
[4,] -493.2138532 -1630.7839588
[5,] -1135.9733763 -493.2138532
[6,] -1342.5736137 -1135.9733763
[7,] -468.3258422 -1342.5736137
[8,] -495.9263473 -468.3258422
[9,] 84.5248815 -495.9263473
[10,] 1455.1956751 84.5248815
[11,] -877.4694952 1455.1956751
[12,] 196.4237011 -877.4694952
[13,] 657.6933635 196.4237011
[14,] 928.9413763 657.6933635
[15,] 100.7743347 928.9413763
[16,] -347.6824620 100.7743347
[17,] -1195.9904081 -347.6824620
[18,] 516.8739375 -1195.9904081
[19,] 877.5651792 516.8739375
[20,] -323.2209084 877.5651792
[21,] 1614.9509795 -323.2209084
[22,] 2348.6238067 1614.9509795
[23,] 1008.5183688 2348.6238067
[24,] 498.5223432 1008.5183688
[25,] 1478.3519601 498.5223432
[26,] -1023.0993050 1478.3519601
[27,] 1913.3722253 -1023.0993050
[28,] 86.6665761 1913.3722253
[29,] -719.5807012 86.6665761
[30,] 315.5728749 -719.5807012
[31,] 257.7689434 315.5728749
[32,] 670.0639648 257.7689434
[33,] 578.2918592 670.0639648
[34,] -468.1818801 578.2918592
[35,] 1247.4117546 -468.1818801
[36,] 399.3430621 1247.4117546
[37,] 1024.8346851 399.3430621
[38,] -1618.7742441 1024.8346851
[39,] 235.8907471 -1618.7742441
[40,] -603.1761451 235.8907471
[41,] 442.9804375 -603.1761451
[42,] 328.0995669 442.9804375
[43,] -872.3461908 328.0995669
[44,] 1318.0171882 -872.3461908
[45,] 844.0728174 1318.0171882
[46,] -1590.6423445 844.0728174
[47,] 1046.1574729 -1590.6423445
[48,] 1007.8906812 1046.1574729
[49,] -1412.7522133 1007.8906812
[50,] 2154.6523418 -1412.7522133
[51,] -288.7727098 2154.6523418
[52,] 1115.6772577 -288.7727098
[53,] 879.1719557 1115.6772577
[54,] -1228.5899052 879.1719557
[55,] 205.3379104 -1228.5899052
[56,] -1168.9338974 205.3379104
[57,] -3121.8405376 -1168.9338974
[58,] -1744.9952573 -3121.8405376
[59,] -2424.6181011 -1744.9952573
[60,] -804.2731005 -2424.6181011
[61,] -1208.9341113 -804.2731005
[62,] 0.7352401 -1208.9341113
[63,] -330.4806386 0.7352401
[64,] 241.7286265 -330.4806386
[65,] 1729.3920925 241.7286265
[66,] 1410.6171396 1729.3920925
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -539.1936842 -1297.9066872
2 -442.4554092 -539.1936842
3 -1630.7839588 -442.4554092
4 -493.2138532 -1630.7839588
5 -1135.9733763 -493.2138532
6 -1342.5736137 -1135.9733763
7 -468.3258422 -1342.5736137
8 -495.9263473 -468.3258422
9 84.5248815 -495.9263473
10 1455.1956751 84.5248815
11 -877.4694952 1455.1956751
12 196.4237011 -877.4694952
13 657.6933635 196.4237011
14 928.9413763 657.6933635
15 100.7743347 928.9413763
16 -347.6824620 100.7743347
17 -1195.9904081 -347.6824620
18 516.8739375 -1195.9904081
19 877.5651792 516.8739375
20 -323.2209084 877.5651792
21 1614.9509795 -323.2209084
22 2348.6238067 1614.9509795
23 1008.5183688 2348.6238067
24 498.5223432 1008.5183688
25 1478.3519601 498.5223432
26 -1023.0993050 1478.3519601
27 1913.3722253 -1023.0993050
28 86.6665761 1913.3722253
29 -719.5807012 86.6665761
30 315.5728749 -719.5807012
31 257.7689434 315.5728749
32 670.0639648 257.7689434
33 578.2918592 670.0639648
34 -468.1818801 578.2918592
35 1247.4117546 -468.1818801
36 399.3430621 1247.4117546
37 1024.8346851 399.3430621
38 -1618.7742441 1024.8346851
39 235.8907471 -1618.7742441
40 -603.1761451 235.8907471
41 442.9804375 -603.1761451
42 328.0995669 442.9804375
43 -872.3461908 328.0995669
44 1318.0171882 -872.3461908
45 844.0728174 1318.0171882
46 -1590.6423445 844.0728174
47 1046.1574729 -1590.6423445
48 1007.8906812 1046.1574729
49 -1412.7522133 1007.8906812
50 2154.6523418 -1412.7522133
51 -288.7727098 2154.6523418
52 1115.6772577 -288.7727098
53 879.1719557 1115.6772577
54 -1228.5899052 879.1719557
55 205.3379104 -1228.5899052
56 -1168.9338974 205.3379104
57 -3121.8405376 -1168.9338974
58 -1744.9952573 -3121.8405376
59 -2424.6181011 -1744.9952573
60 -804.2731005 -2424.6181011
61 -1208.9341113 -804.2731005
62 0.7352401 -1208.9341113
63 -330.4806386 0.7352401
64 241.7286265 -330.4806386
65 1729.3920925 241.7286265
66 1410.6171396 1729.3920925
> 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/7epo61258477215.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/84dn41258477215.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/9yopk1258477215.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/10u6141258477215.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/11f8v11258477215.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/12qpr21258477215.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/13sxur1258477215.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/145rpm1258477215.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/15tlcm1258477215.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/16rna71258477215.tab")
+ }
>
> system("convert tmp/1ty3q1258477215.ps tmp/1ty3q1258477215.png")
> system("convert tmp/2sf9t1258477215.ps tmp/2sf9t1258477215.png")
> system("convert tmp/39qno1258477215.ps tmp/39qno1258477215.png")
> system("convert tmp/4c5mh1258477215.ps tmp/4c5mh1258477215.png")
> system("convert tmp/5n9r41258477215.ps tmp/5n9r41258477215.png")
> system("convert tmp/6v3f11258477215.ps tmp/6v3f11258477215.png")
> system("convert tmp/7epo61258477215.ps tmp/7epo61258477215.png")
> system("convert tmp/84dn41258477215.ps tmp/84dn41258477215.png")
> system("convert tmp/9yopk1258477215.ps tmp/9yopk1258477215.png")
> system("convert tmp/10u6141258477215.ps tmp/10u6141258477215.png")
>
>
> proc.time()
user system elapsed
2.573 1.622 4.168