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
+ ,22686.8
+ ,2.01
+ ,19435.1
+ ,18714.9
+ ,20396.7
+ ,2.01
+ ,22686.8
+ ,19435.1
+ ,19233.6
+ ,2.01
+ ,20396.7
+ ,22686.8
+ ,22751
+ ,2.01
+ ,19233.6
+ ,20396.7
+ ,19864
+ ,2.01
+ ,22751
+ ,19233.6
+ ,17165.4
+ ,2.02
+ ,19864
+ ,22751
+ ,22309.7
+ ,2.02
+ ,17165.4
+ ,19864
+ ,21786.3
+ ,2.03
+ ,22309.7
+ ,17165.4
+ ,21927.6
+ ,2.05
+ ,21786.3
+ ,22309.7
+ ,20957.9
+ ,2.08
+ ,21927.6
+ ,21786.3
+ ,19726
+ ,2.07
+ ,20957.9
+ ,21927.6
+ ,21315.7
+ ,2.06
+ ,19726
+ ,20957.9
+ ,24771.5
+ ,2.05
+ ,21315.7
+ ,19726
+ ,22592.4
+ ,2.05
+ ,24771.5
+ ,21315.7
+ ,21942.1
+ ,2.05
+ ,22592.4
+ ,24771.5
+ ,23973.7
+ ,2.05
+ ,21942.1
+ ,22592.4
+ ,20815.7
+ ,2.05
+ ,23973.7
+ ,21942.1
+ ,19931.4
+ ,2.06
+ ,20815.7
+ ,23973.7
+ ,24436.8
+ ,2.06
+ ,19931.4
+ ,20815.7
+ ,22838.7
+ ,2.07
+ ,24436.8
+ ,19931.4
+ ,24465.3
+ ,2.07
+ ,22838.7
+ ,24436.8
+ ,23007.3
+ ,2.3
+ ,24465.3
+ ,22838.7
+ ,22720.8
+ ,2.31
+ ,23007.3
+ ,24465.3
+ ,23045.7
+ ,2.31
+ ,22720.8
+ ,23007.3
+ ,27198.5
+ ,2.53
+ ,23045.7
+ ,22720.8
+ ,22401.9
+ ,2.58
+ ,27198.5
+ ,23045.7
+ ,25122.7
+ ,2.59
+ ,22401.9
+ ,27198.5
+ ,26100.5
+ ,2.73
+ ,25122.7
+ ,22401.9
+ ,22904.9
+ ,2.82
+ ,26100.5
+ ,25122.7
+ ,22040.4
+ ,3
+ ,22904.9
+ ,26100.5
+ ,25981.5
+ ,3.04
+ ,22040.4
+ ,22904.9
+ ,26157.1
+ ,3.23
+ ,25981.5
+ ,22040.4
+ ,25975.4
+ ,3.32
+ ,26157.1
+ ,25981.5
+ ,22589.8
+ ,3.49
+ ,25975.4
+ ,26157.1
+ ,25370.4
+ ,3.57
+ ,22589.8
+ ,25975.4
+ ,25091.1
+ ,3.56
+ ,25370.4
+ ,22589.8
+ ,28760.9
+ ,3.72
+ ,25091.1
+ ,25370.4
+ ,24325.9
+ ,3.82
+ ,28760.9
+ ,25091.1
+ ,25821.7
+ ,3.82
+ ,24325.9
+ ,28760.9
+ ,27645.7
+ ,3.98
+ ,25821.7
+ ,24325.9
+ ,26296.9
+ ,4.06
+ ,27645.7
+ ,25821.7
+ ,24141.5
+ ,4.08
+ ,26296.9
+ ,27645.7
+ ,27268.1
+ ,4.19
+ ,24141.5
+ ,26296.9
+ ,29060.3
+ ,4.16
+ ,27268.1
+ ,24141.5
+ ,28226.4
+ ,4.17
+ ,29060.3
+ ,27268.1
+ ,23268.5
+ ,4.21
+ ,28226.4
+ ,29060.3
+ ,26938.2
+ ,4.21
+ ,23268.5
+ ,28226.4
+ ,27217.5
+ ,4.17
+ ,26938.2
+ ,23268.5
+ ,27540.5
+ ,4.19
+ ,27217.5
+ ,26938.2
+ ,29167.6
+ ,4.25
+ ,27540.5
+ ,27217.5
+ ,26671.5
+ ,4.25
+ ,29167.6
+ ,27540.5
+ ,30184
+ ,4.2
+ ,26671.5
+ ,29167.6
+ ,28422.3
+ ,4.33
+ ,30184
+ ,26671.5
+ ,23774.3
+ ,4.41
+ ,28422.3
+ ,30184
+ ,29601
+ ,4.56
+ ,23774.3
+ ,28422.3
+ ,28523.6
+ ,5.18
+ ,29601
+ ,23774.3
+ ,23622
+ ,3.42
+ ,28523.6
+ ,29601
+ ,21320.3
+ ,2.71
+ ,23622
+ ,28523.6
+ ,20423.6
+ ,2.29
+ ,21320.3
+ ,23622
+ ,21174.9
+ ,2
+ ,20423.6
+ ,21320.3
+ ,23050.2
+ ,1.64
+ ,21174.9
+ ,20423.6
+ ,21202.9
+ ,1.3
+ ,23050.2
+ ,21174.9
+ ,20476.4
+ ,1.08
+ ,21202.9
+ ,23050.2
+ ,23173.3
+ ,1
+ ,20476.4
+ ,21202.9
+ ,22468
+ ,1
+ ,23173.3
+ ,20476.4
+ ,19842.7
+ ,1
+ ,22468
+ ,23173.3)
+ ,dim=c(4
+ ,67)
+ ,dimnames=list(c('Y'
+ ,'X'
+ ,'Y1'
+ ,'Y2')
+ ,1:67))
> y <- array(NA,dim=c(4,67),dimnames=list(c('Y','X','Y1','Y2'),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 M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 19435.1 2.01 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 0 1 0 0 0 0 0 0 0 0 0 2
3 20396.7 2.01 22686.8 19435.1 0 0 1 0 0 0 0 0 0 0 0 3
4 19233.6 2.01 20396.7 22686.8 0 0 0 1 0 0 0 0 0 0 0 4
5 22751.0 2.01 19233.6 20396.7 0 0 0 0 1 0 0 0 0 0 0 5
6 19864.0 2.01 22751.0 19233.6 0 0 0 0 0 1 0 0 0 0 0 6
7 17165.4 2.02 19864.0 22751.0 0 0 0 0 0 0 1 0 0 0 0 7
8 22309.7 2.02 17165.4 19864.0 0 0 0 0 0 0 0 1 0 0 0 8
9 21786.3 2.03 22309.7 17165.4 0 0 0 0 0 0 0 0 1 0 0 9
10 21927.6 2.05 21786.3 22309.7 0 0 0 0 0 0 0 0 0 1 0 10
11 20957.9 2.08 21927.6 21786.3 0 0 0 0 0 0 0 0 0 0 1 11
12 19726.0 2.07 20957.9 21927.6 0 0 0 0 0 0 0 0 0 0 0 12
13 21315.7 2.06 19726.0 20957.9 1 0 0 0 0 0 0 0 0 0 0 13
14 24771.5 2.05 21315.7 19726.0 0 1 0 0 0 0 0 0 0 0 0 14
15 22592.4 2.05 24771.5 21315.7 0 0 1 0 0 0 0 0 0 0 0 15
16 21942.1 2.05 22592.4 24771.5 0 0 0 1 0 0 0 0 0 0 0 16
17 23973.7 2.05 21942.1 22592.4 0 0 0 0 1 0 0 0 0 0 0 17
18 20815.7 2.05 23973.7 21942.1 0 0 0 0 0 1 0 0 0 0 0 18
19 19931.4 2.06 20815.7 23973.7 0 0 0 0 0 0 1 0 0 0 0 19
20 24436.8 2.06 19931.4 20815.7 0 0 0 0 0 0 0 1 0 0 0 20
21 22838.7 2.07 24436.8 19931.4 0 0 0 0 0 0 0 0 1 0 0 21
22 24465.3 2.07 22838.7 24436.8 0 0 0 0 0 0 0 0 0 1 0 22
23 23007.3 2.30 24465.3 22838.7 0 0 0 0 0 0 0 0 0 0 1 23
24 22720.8 2.31 23007.3 24465.3 0 0 0 0 0 0 0 0 0 0 0 24
25 23045.7 2.31 22720.8 23007.3 1 0 0 0 0 0 0 0 0 0 0 25
26 27198.5 2.53 23045.7 22720.8 0 1 0 0 0 0 0 0 0 0 0 26
27 22401.9 2.58 27198.5 23045.7 0 0 1 0 0 0 0 0 0 0 0 27
28 25122.7 2.59 22401.9 27198.5 0 0 0 1 0 0 0 0 0 0 0 28
29 26100.5 2.73 25122.7 22401.9 0 0 0 0 1 0 0 0 0 0 0 29
30 22904.9 2.82 26100.5 25122.7 0 0 0 0 0 1 0 0 0 0 0 30
31 22040.4 3.00 22904.9 26100.5 0 0 0 0 0 0 1 0 0 0 0 31
32 25981.5 3.04 22040.4 22904.9 0 0 0 0 0 0 0 1 0 0 0 32
33 26157.1 3.23 25981.5 22040.4 0 0 0 0 0 0 0 0 1 0 0 33
34 25975.4 3.32 26157.1 25981.5 0 0 0 0 0 0 0 0 0 1 0 34
35 22589.8 3.49 25975.4 26157.1 0 0 0 0 0 0 0 0 0 0 1 35
36 25370.4 3.57 22589.8 25975.4 0 0 0 0 0 0 0 0 0 0 0 36
37 25091.1 3.56 25370.4 22589.8 1 0 0 0 0 0 0 0 0 0 0 37
38 28760.9 3.72 25091.1 25370.4 0 1 0 0 0 0 0 0 0 0 0 38
39 24325.9 3.82 28760.9 25091.1 0 0 1 0 0 0 0 0 0 0 0 39
40 25821.7 3.82 24325.9 28760.9 0 0 0 1 0 0 0 0 0 0 0 40
41 27645.7 3.98 25821.7 24325.9 0 0 0 0 1 0 0 0 0 0 0 41
42 26296.9 4.06 27645.7 25821.7 0 0 0 0 0 1 0 0 0 0 0 42
43 24141.5 4.08 26296.9 27645.7 0 0 0 0 0 0 1 0 0 0 0 43
44 27268.1 4.19 24141.5 26296.9 0 0 0 0 0 0 0 1 0 0 0 44
45 29060.3 4.16 27268.1 24141.5 0 0 0 0 0 0 0 0 1 0 0 45
46 28226.4 4.17 29060.3 27268.1 0 0 0 0 0 0 0 0 0 1 0 46
47 23268.5 4.21 28226.4 29060.3 0 0 0 0 0 0 0 0 0 0 1 47
48 26938.2 4.21 23268.5 28226.4 0 0 0 0 0 0 0 0 0 0 0 48
49 27217.5 4.17 26938.2 23268.5 1 0 0 0 0 0 0 0 0 0 0 49
50 27540.5 4.19 27217.5 26938.2 0 1 0 0 0 0 0 0 0 0 0 50
51 29167.6 4.25 27540.5 27217.5 0 0 1 0 0 0 0 0 0 0 0 51
52 26671.5 4.25 29167.6 27540.5 0 0 0 1 0 0 0 0 0 0 0 52
53 30184.0 4.20 26671.5 29167.6 0 0 0 0 1 0 0 0 0 0 0 53
54 28422.3 4.33 30184.0 26671.5 0 0 0 0 0 1 0 0 0 0 0 54
55 23774.3 4.41 28422.3 30184.0 0 0 0 0 0 0 1 0 0 0 0 55
56 29601.0 4.56 23774.3 28422.3 0 0 0 0 0 0 0 1 0 0 0 56
57 28523.6 5.18 29601.0 23774.3 0 0 0 0 0 0 0 0 1 0 0 57
58 23622.0 3.42 28523.6 29601.0 0 0 0 0 0 0 0 0 0 1 0 58
59 21320.3 2.71 23622.0 28523.6 0 0 0 0 0 0 0 0 0 0 1 59
60 20423.6 2.29 21320.3 23622.0 0 0 0 0 0 0 0 0 0 0 0 60
61 21174.9 2.00 20423.6 21320.3 1 0 0 0 0 0 0 0 0 0 0 61
62 23050.2 1.64 21174.9 20423.6 0 1 0 0 0 0 0 0 0 0 0 62
63 21202.9 1.30 23050.2 21174.9 0 0 1 0 0 0 0 0 0 0 0 63
64 20476.4 1.08 21202.9 23050.2 0 0 0 1 0 0 0 0 0 0 0 64
65 23173.3 1.00 20476.4 21202.9 0 0 0 0 1 0 0 0 0 0 0 65
66 22468.0 1.00 23173.3 20476.4 0 0 0 0 0 1 0 0 0 0 0 66
67 19842.7 1.00 22468.0 23173.3 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 M1 M2
12447.5735 1520.2342 0.1063 0.1257 574.2009 3263.5100
M3 M4 M5 M6 M7 M8
589.2038 382.4990 3052.4960 539.8414 -1942.1129 2766.2863
M9 M10 M11 t
2060.4529 1169.1283 -1273.0040 19.7034
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-3089.3 -916.5 125.4 810.4 2411.7
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 12447.5735 3567.0011 3.490 0.001006 **
X 1520.2342 445.6710 3.411 0.001274 **
Y1 0.1063 0.1556 0.683 0.497446
Y2 0.1257 0.1487 0.845 0.402073
M1 574.2009 880.3050 0.652 0.517153
M2 3263.5100 869.6283 3.753 0.000449 ***
M3 589.2038 1034.2543 0.570 0.571389
M4 382.4990 859.6450 0.445 0.658239
M5 3052.4960 841.8480 3.626 0.000665 ***
M6 539.8414 1009.8230 0.535 0.595256
M7 -1942.1129 837.0093 -2.320 0.024362 *
M8 2766.2863 883.1654 3.132 0.002872 **
M9 2060.4529 1143.2413 1.802 0.077411 .
M10 1169.1283 980.1189 1.193 0.238451
M11 -1273.0040 920.7473 -1.383 0.172819
t 19.7034 11.7311 1.680 0.099154 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1310 on 51 degrees of freedom
Multiple R-squared: 0.8525, Adjusted R-squared: 0.8091
F-statistic: 19.65 on 15 and 51 DF, p-value: 4.207e-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.0846041860 0.1692083720 0.9153958
[2,] 0.0307778020 0.0615556040 0.9692222
[3,] 0.0103109237 0.0206218474 0.9896891
[4,] 0.0042656070 0.0085312141 0.9957344
[5,] 0.0057599020 0.0115198041 0.9942401
[6,] 0.0057648283 0.0115296566 0.9942352
[7,] 0.0036566175 0.0073132349 0.9963434
[8,] 0.0018313112 0.0036626225 0.9981687
[9,] 0.0133501745 0.0267003491 0.9866498
[10,] 0.0168706877 0.0337413755 0.9831293
[11,] 0.0095546003 0.0191092006 0.9904454
[12,] 0.0080533211 0.0161066421 0.9919467
[13,] 0.0044247754 0.0088495507 0.9955752
[14,] 0.0020640294 0.0041280588 0.9979360
[15,] 0.0010899808 0.0021799615 0.9989100
[16,] 0.0004743490 0.0009486981 0.9995257
[17,] 0.0030201021 0.0060402041 0.9969799
[18,] 0.0017985481 0.0035970962 0.9982015
[19,] 0.0008098961 0.0016197921 0.9991901
[20,] 0.0006798304 0.0013596609 0.9993202
[21,] 0.0012257110 0.0024514219 0.9987743
[22,] 0.0005651463 0.0011302925 0.9994349
[23,] 0.0004164796 0.0008329592 0.9995835
[24,] 0.0019736931 0.0039473861 0.9980263
[25,] 0.0060478313 0.0120956626 0.9939522
[26,] 0.0252664069 0.0505328139 0.9747336
[27,] 0.0227075486 0.0454150972 0.9772925
[28,] 0.0118362353 0.0236724707 0.9881638
[29,] 0.1267450791 0.2534901581 0.8732549
[30,] 0.0702169231 0.1404338462 0.9297831
> postscript(file="/var/www/html/rcomp/tmp/1n5q51258477094.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/2npd71258477094.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/3l0li1258477094.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/45igl1258477094.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/56gou1258477094.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 6
-1442.13487 -537.63329 -609.38780 -1750.58657 -511.43902 -1133.34169
7 8 9 10 11 12
-1519.91214 -453.98886 -514.35400 -122.61714 1335.25054 -1088.80102
13 14 15 16 17 18
175.03884 922.79476 831.07660 165.22784 -149.89947 -949.25389
19 20 21 22 23 24
693.99718 962.15783 -332.95751 1769.04232 2411.67649 767.89876
25 26 27 28 29 30
712.57068 1823.36124 -877.04650 2003.73897 392.44142 -892.89516
31 32 33 34 35 36
648.13477 293.80551 556.26259 595.45185 -628.90473 1120.19616
37 38 39 40 41 42
391.96472 789.80246 -1497.73523 195.49710 -515.18823 125.43379
43 44 45 46 47 48
316.09532 -1053.95661 1408.37332 847.43662 -1885.38250 1123.57705
49 50 51 52 53 54
1102.58905 -1804.66138 2316.38550 -206.31271 753.45013 1227.24296
55 56 57 58 59 60
-1334.18255 251.98213 -1117.32440 -3089.31366 -1232.63981 -1922.87094
61 62 63 64 65 66
-940.02842 -1193.66379 -163.29256 -407.56464 30.63519 1622.81399
67
1195.86742
> postscript(file="/var/www/html/rcomp/tmp/6b7b71258477094.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 -1442.13487 NA
1 -537.63329 -1442.13487
2 -609.38780 -537.63329
3 -1750.58657 -609.38780
4 -511.43902 -1750.58657
5 -1133.34169 -511.43902
6 -1519.91214 -1133.34169
7 -453.98886 -1519.91214
8 -514.35400 -453.98886
9 -122.61714 -514.35400
10 1335.25054 -122.61714
11 -1088.80102 1335.25054
12 175.03884 -1088.80102
13 922.79476 175.03884
14 831.07660 922.79476
15 165.22784 831.07660
16 -149.89947 165.22784
17 -949.25389 -149.89947
18 693.99718 -949.25389
19 962.15783 693.99718
20 -332.95751 962.15783
21 1769.04232 -332.95751
22 2411.67649 1769.04232
23 767.89876 2411.67649
24 712.57068 767.89876
25 1823.36124 712.57068
26 -877.04650 1823.36124
27 2003.73897 -877.04650
28 392.44142 2003.73897
29 -892.89516 392.44142
30 648.13477 -892.89516
31 293.80551 648.13477
32 556.26259 293.80551
33 595.45185 556.26259
34 -628.90473 595.45185
35 1120.19616 -628.90473
36 391.96472 1120.19616
37 789.80246 391.96472
38 -1497.73523 789.80246
39 195.49710 -1497.73523
40 -515.18823 195.49710
41 125.43379 -515.18823
42 316.09532 125.43379
43 -1053.95661 316.09532
44 1408.37332 -1053.95661
45 847.43662 1408.37332
46 -1885.38250 847.43662
47 1123.57705 -1885.38250
48 1102.58905 1123.57705
49 -1804.66138 1102.58905
50 2316.38550 -1804.66138
51 -206.31271 2316.38550
52 753.45013 -206.31271
53 1227.24296 753.45013
54 -1334.18255 1227.24296
55 251.98213 -1334.18255
56 -1117.32440 251.98213
57 -3089.31366 -1117.32440
58 -1232.63981 -3089.31366
59 -1922.87094 -1232.63981
60 -940.02842 -1922.87094
61 -1193.66379 -940.02842
62 -163.29256 -1193.66379
63 -407.56464 -163.29256
64 30.63519 -407.56464
65 1622.81399 30.63519
66 1195.86742 1622.81399
67 NA 1195.86742
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -537.63329 -1442.13487
[2,] -609.38780 -537.63329
[3,] -1750.58657 -609.38780
[4,] -511.43902 -1750.58657
[5,] -1133.34169 -511.43902
[6,] -1519.91214 -1133.34169
[7,] -453.98886 -1519.91214
[8,] -514.35400 -453.98886
[9,] -122.61714 -514.35400
[10,] 1335.25054 -122.61714
[11,] -1088.80102 1335.25054
[12,] 175.03884 -1088.80102
[13,] 922.79476 175.03884
[14,] 831.07660 922.79476
[15,] 165.22784 831.07660
[16,] -149.89947 165.22784
[17,] -949.25389 -149.89947
[18,] 693.99718 -949.25389
[19,] 962.15783 693.99718
[20,] -332.95751 962.15783
[21,] 1769.04232 -332.95751
[22,] 2411.67649 1769.04232
[23,] 767.89876 2411.67649
[24,] 712.57068 767.89876
[25,] 1823.36124 712.57068
[26,] -877.04650 1823.36124
[27,] 2003.73897 -877.04650
[28,] 392.44142 2003.73897
[29,] -892.89516 392.44142
[30,] 648.13477 -892.89516
[31,] 293.80551 648.13477
[32,] 556.26259 293.80551
[33,] 595.45185 556.26259
[34,] -628.90473 595.45185
[35,] 1120.19616 -628.90473
[36,] 391.96472 1120.19616
[37,] 789.80246 391.96472
[38,] -1497.73523 789.80246
[39,] 195.49710 -1497.73523
[40,] -515.18823 195.49710
[41,] 125.43379 -515.18823
[42,] 316.09532 125.43379
[43,] -1053.95661 316.09532
[44,] 1408.37332 -1053.95661
[45,] 847.43662 1408.37332
[46,] -1885.38250 847.43662
[47,] 1123.57705 -1885.38250
[48,] 1102.58905 1123.57705
[49,] -1804.66138 1102.58905
[50,] 2316.38550 -1804.66138
[51,] -206.31271 2316.38550
[52,] 753.45013 -206.31271
[53,] 1227.24296 753.45013
[54,] -1334.18255 1227.24296
[55,] 251.98213 -1334.18255
[56,] -1117.32440 251.98213
[57,] -3089.31366 -1117.32440
[58,] -1232.63981 -3089.31366
[59,] -1922.87094 -1232.63981
[60,] -940.02842 -1922.87094
[61,] -1193.66379 -940.02842
[62,] -163.29256 -1193.66379
[63,] -407.56464 -163.29256
[64,] 30.63519 -407.56464
[65,] 1622.81399 30.63519
[66,] 1195.86742 1622.81399
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -537.63329 -1442.13487
2 -609.38780 -537.63329
3 -1750.58657 -609.38780
4 -511.43902 -1750.58657
5 -1133.34169 -511.43902
6 -1519.91214 -1133.34169
7 -453.98886 -1519.91214
8 -514.35400 -453.98886
9 -122.61714 -514.35400
10 1335.25054 -122.61714
11 -1088.80102 1335.25054
12 175.03884 -1088.80102
13 922.79476 175.03884
14 831.07660 922.79476
15 165.22784 831.07660
16 -149.89947 165.22784
17 -949.25389 -149.89947
18 693.99718 -949.25389
19 962.15783 693.99718
20 -332.95751 962.15783
21 1769.04232 -332.95751
22 2411.67649 1769.04232
23 767.89876 2411.67649
24 712.57068 767.89876
25 1823.36124 712.57068
26 -877.04650 1823.36124
27 2003.73897 -877.04650
28 392.44142 2003.73897
29 -892.89516 392.44142
30 648.13477 -892.89516
31 293.80551 648.13477
32 556.26259 293.80551
33 595.45185 556.26259
34 -628.90473 595.45185
35 1120.19616 -628.90473
36 391.96472 1120.19616
37 789.80246 391.96472
38 -1497.73523 789.80246
39 195.49710 -1497.73523
40 -515.18823 195.49710
41 125.43379 -515.18823
42 316.09532 125.43379
43 -1053.95661 316.09532
44 1408.37332 -1053.95661
45 847.43662 1408.37332
46 -1885.38250 847.43662
47 1123.57705 -1885.38250
48 1102.58905 1123.57705
49 -1804.66138 1102.58905
50 2316.38550 -1804.66138
51 -206.31271 2316.38550
52 753.45013 -206.31271
53 1227.24296 753.45013
54 -1334.18255 1227.24296
55 251.98213 -1334.18255
56 -1117.32440 251.98213
57 -3089.31366 -1117.32440
58 -1232.63981 -3089.31366
59 -1922.87094 -1232.63981
60 -940.02842 -1922.87094
61 -1193.66379 -940.02842
62 -163.29256 -1193.66379
63 -407.56464 -163.29256
64 30.63519 -407.56464
65 1622.81399 30.63519
66 1195.86742 1622.81399
> 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/7hany1258477094.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/8km9i1258477094.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/9c8df1258477094.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/10hg4f1258477094.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/11xf8r1258477094.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/12n2cs1258477094.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/1390ea1258477094.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/14rvdd1258477095.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/151lrk1258477095.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/16ehft1258477095.tab")
+ }
>
> system("convert tmp/1n5q51258477094.ps tmp/1n5q51258477094.png")
> system("convert tmp/2npd71258477094.ps tmp/2npd71258477094.png")
> system("convert tmp/3l0li1258477094.ps tmp/3l0li1258477094.png")
> system("convert tmp/45igl1258477094.ps tmp/45igl1258477094.png")
> system("convert tmp/56gou1258477094.ps tmp/56gou1258477094.png")
> system("convert tmp/6b7b71258477094.ps tmp/6b7b71258477094.png")
> system("convert tmp/7hany1258477094.ps tmp/7hany1258477094.png")
> system("convert tmp/8km9i1258477094.ps tmp/8km9i1258477094.png")
> system("convert tmp/9c8df1258477094.ps tmp/9c8df1258477094.png")
> system("convert tmp/10hg4f1258477094.ps tmp/10hg4f1258477094.png")
>
>
> proc.time()
user system elapsed
2.629 1.635 3.829