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(2.6,30.5,2.4,28.6,2.5,30,2.7,28.2,3.2,27.6,2.8,24.9,2.8,23.8,3,24.3,3.1,23.6,3.1,24.2,3,28.1,2.4,30.1,2.7,31.1,3,32,2.7,32.4,2.7,34,2,35.1,2.4,37.1,2.6,37.3,2.4,38.1,2.3,39.5,2.4,38.3,2.5,37.3,2.6,38.7,2.6,37.5,2.6,38.7,2.7,37.9,2.8,36.6,2.6,35.5,2.6,37.6,2,38.6,2,40.3,2.1,39,1.9,36.8,2,36.5,2.5,34.1,2.9,34.2,3.3,31.9,3.5,33.7,3.8,33.5,4.6,33.8,4.4,29.9,5.3,32.3,5.8,30.5,5.9,28.5,5.6,29,5.8,23.8,5.5,17.9,4.6,9.9,4.2,3,4,4.2,3.5,0.4,2.3,0,2.2,2.4,1.4,4.2,0.6,8.2,0,9,0.5,13.6,0.1,14,0.1,17.6),dim=c(2,60),dimnames=list(c('Y','X'),1:60))
> y <- array(NA,dim=c(2,60),dimnames=list(c('Y','X'),1:60))
> for (i in 1:dim(x)[1])
+ {
+ for (j in 1:dim(x)[2])
+ {
+ y[i,j] <- as.numeric(x[i,j])
+ }
+ }
> par3 = 'Linear Trend'
> par2 = 'Include Monthly Dummies'
> par1 = '1'
> #'GNU S' R Code compiled by R2WASP v. 1.0.44 ()
> #Author: Prof. Dr. P. Wessa
> #To cite this work: AUTHOR(S), (YEAR), YOUR SOFTWARE TITLE (vNUMBER) in Free Statistics Software (v$_version), Office for Research Development and Education, URL http://www.wessa.net/rwasp_YOURPAGE.wasp/
> #Source of accompanying publication: Office for Research, Development, and Education
> #Technical description: Write here your technical program description (don't use hard returns!)
> library(lattice)
> library(lmtest)
Loading required package: zoo
Attaching package: 'zoo'
The following object(s) are masked from package:base :
as.Date.numeric
> n25 <- 25 #minimum number of obs. for Goldfeld-Quandt test
> par1 <- as.numeric(par1)
> x <- t(y)
> k <- length(x[1,])
> n <- length(x[,1])
> x1 <- cbind(x[,par1], x[,1:k!=par1])
> mycolnames <- c(colnames(x)[par1], colnames(x)[1:k!=par1])
> colnames(x1) <- mycolnames #colnames(x)[par1]
> x <- x1
> if (par3 == 'First Differences'){
+ x2 <- array(0, dim=c(n-1,k), dimnames=list(1:(n-1), paste('(1-B)',colnames(x),sep='')))
+ for (i in 1:n-1) {
+ for (j in 1:k) {
+ x2[i,j] <- x[i+1,j] - x[i,j]
+ }
+ }
+ x <- x2
+ }
> if (par2 == 'Include Monthly Dummies'){
+ x2 <- array(0, dim=c(n,11), dimnames=list(1:n, paste('M', seq(1:11), sep ='')))
+ for (i in 1:11){
+ x2[seq(i,n,12),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> if (par2 == 'Include Quarterly Dummies'){
+ x2 <- array(0, dim=c(n,3), dimnames=list(1:n, paste('Q', seq(1:3), sep ='')))
+ for (i in 1:3){
+ x2[seq(i,n,4),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> k <- length(x[1,])
> if (par3 == 'Linear Trend'){
+ x <- cbind(x, c(1:n))
+ colnames(x)[k+1] <- 't'
+ }
> x
Y X M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 2.6 30.5 1 0 0 0 0 0 0 0 0 0 0 1
2 2.4 28.6 0 1 0 0 0 0 0 0 0 0 0 2
3 2.5 30.0 0 0 1 0 0 0 0 0 0 0 0 3
4 2.7 28.2 0 0 0 1 0 0 0 0 0 0 0 4
5 3.2 27.6 0 0 0 0 1 0 0 0 0 0 0 5
6 2.8 24.9 0 0 0 0 0 1 0 0 0 0 0 6
7 2.8 23.8 0 0 0 0 0 0 1 0 0 0 0 7
8 3.0 24.3 0 0 0 0 0 0 0 1 0 0 0 8
9 3.1 23.6 0 0 0 0 0 0 0 0 1 0 0 9
10 3.1 24.2 0 0 0 0 0 0 0 0 0 1 0 10
11 3.0 28.1 0 0 0 0 0 0 0 0 0 0 1 11
12 2.4 30.1 0 0 0 0 0 0 0 0 0 0 0 12
13 2.7 31.1 1 0 0 0 0 0 0 0 0 0 0 13
14 3.0 32.0 0 1 0 0 0 0 0 0 0 0 0 14
15 2.7 32.4 0 0 1 0 0 0 0 0 0 0 0 15
16 2.7 34.0 0 0 0 1 0 0 0 0 0 0 0 16
17 2.0 35.1 0 0 0 0 1 0 0 0 0 0 0 17
18 2.4 37.1 0 0 0 0 0 1 0 0 0 0 0 18
19 2.6 37.3 0 0 0 0 0 0 1 0 0 0 0 19
20 2.4 38.1 0 0 0 0 0 0 0 1 0 0 0 20
21 2.3 39.5 0 0 0 0 0 0 0 0 1 0 0 21
22 2.4 38.3 0 0 0 0 0 0 0 0 0 1 0 22
23 2.5 37.3 0 0 0 0 0 0 0 0 0 0 1 23
24 2.6 38.7 0 0 0 0 0 0 0 0 0 0 0 24
25 2.6 37.5 1 0 0 0 0 0 0 0 0 0 0 25
26 2.6 38.7 0 1 0 0 0 0 0 0 0 0 0 26
27 2.7 37.9 0 0 1 0 0 0 0 0 0 0 0 27
28 2.8 36.6 0 0 0 1 0 0 0 0 0 0 0 28
29 2.6 35.5 0 0 0 0 1 0 0 0 0 0 0 29
30 2.6 37.6 0 0 0 0 0 1 0 0 0 0 0 30
31 2.0 38.6 0 0 0 0 0 0 1 0 0 0 0 31
32 2.0 40.3 0 0 0 0 0 0 0 1 0 0 0 32
33 2.1 39.0 0 0 0 0 0 0 0 0 1 0 0 33
34 1.9 36.8 0 0 0 0 0 0 0 0 0 1 0 34
35 2.0 36.5 0 0 0 0 0 0 0 0 0 0 1 35
36 2.5 34.1 0 0 0 0 0 0 0 0 0 0 0 36
37 2.9 34.2 1 0 0 0 0 0 0 0 0 0 0 37
38 3.3 31.9 0 1 0 0 0 0 0 0 0 0 0 38
39 3.5 33.7 0 0 1 0 0 0 0 0 0 0 0 39
40 3.8 33.5 0 0 0 1 0 0 0 0 0 0 0 40
41 4.6 33.8 0 0 0 0 1 0 0 0 0 0 0 41
42 4.4 29.9 0 0 0 0 0 1 0 0 0 0 0 42
43 5.3 32.3 0 0 0 0 0 0 1 0 0 0 0 43
44 5.8 30.5 0 0 0 0 0 0 0 1 0 0 0 44
45 5.9 28.5 0 0 0 0 0 0 0 0 1 0 0 45
46 5.6 29.0 0 0 0 0 0 0 0 0 0 1 0 46
47 5.8 23.8 0 0 0 0 0 0 0 0 0 0 1 47
48 5.5 17.9 0 0 0 0 0 0 0 0 0 0 0 48
49 4.6 9.9 1 0 0 0 0 0 0 0 0 0 0 49
50 4.2 3.0 0 1 0 0 0 0 0 0 0 0 0 50
51 4.0 4.2 0 0 1 0 0 0 0 0 0 0 0 51
52 3.5 0.4 0 0 0 1 0 0 0 0 0 0 0 52
53 2.3 0.0 0 0 0 0 1 0 0 0 0 0 0 53
54 2.2 2.4 0 0 0 0 0 1 0 0 0 0 0 54
55 1.4 4.2 0 0 0 0 0 0 1 0 0 0 0 55
56 0.6 8.2 0 0 0 0 0 0 0 1 0 0 0 56
57 0.0 9.0 0 0 0 0 0 0 0 0 1 0 0 57
58 0.5 13.6 0 0 0 0 0 0 0 0 0 1 0 58
59 0.1 14.0 0 0 0 0 0 0 0 0 0 0 1 59
60 0.1 17.6 0 0 0 0 0 0 0 0 0 0 0 60
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X M1 M2 M3 M4
1.72690 0.01803 0.56307 0.60458 0.55921 0.58810
M5 M6 M7 M8 M9 M10
0.41969 0.34910 0.26265 0.17295 0.08850 0.08926
M11 t
0.06625 0.01094
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-2.6015 -0.5477 -0.2816 0.5757 3.0782
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.72690 1.14073 1.514 0.137
X 0.01803 0.02024 0.891 0.378
M1 0.56307 0.94548 0.596 0.554
M2 0.60458 0.94649 0.639 0.526
M3 0.55921 0.94315 0.593 0.556
M4 0.58810 0.94320 0.624 0.536
M5 0.41969 0.94177 0.446 0.658
M6 0.34910 0.94031 0.371 0.712
M7 0.26265 0.93797 0.280 0.781
M8 0.17295 0.93635 0.185 0.854
M9 0.08850 0.93587 0.095 0.925
M10 0.08926 0.93536 0.095 0.924
M11 0.06625 0.93516 0.071 0.944
t 0.01094 0.01362 0.803 0.426
Residual standard error: 1.478 on 46 degrees of freedom
Multiple R-squared: 0.03745, Adjusted R-squared: -0.2346
F-statistic: 0.1377 on 13 and 46 DF, p-value: 0.9998
> if (n > n25) {
+ kp3 <- k + 3
+ nmkm3 <- n - k - 3
+ gqarr <- array(NA, dim=c(nmkm3-kp3+1,3))
+ numgqtests <- 0
+ numsignificant1 <- 0
+ numsignificant5 <- 0
+ numsignificant10 <- 0
+ for (mypoint in kp3:nmkm3) {
+ j <- 0
+ numgqtests <- numgqtests + 1
+ for (myalt in c('greater', 'two.sided', 'less')) {
+ j <- j + 1
+ gqarr[mypoint-kp3+1,j] <- gqtest(mylm, point=mypoint, alternative=myalt)$p.value
+ }
+ if (gqarr[mypoint-kp3+1,2] < 0.01) numsignificant1 <- numsignificant1 + 1
+ if (gqarr[mypoint-kp3+1,2] < 0.05) numsignificant5 <- numsignificant5 + 1
+ if (gqarr[mypoint-kp3+1,2] < 0.10) numsignificant10 <- numsignificant10 + 1
+ }
+ gqarr
+ }
[,1] [,2] [,3]
[1,] 1.380754e-02 2.761508e-02 0.9861925
[2,] 4.438583e-03 8.877165e-03 0.9955614
[3,] 1.142209e-03 2.284418e-03 0.9988578
[4,] 1.998882e-04 3.997764e-04 0.9998001
[5,] 3.251005e-05 6.502009e-05 0.9999675
[6,] 4.864980e-06 9.729960e-06 0.9999951
[7,] 7.273028e-07 1.454606e-06 0.9999993
[8,] 1.838688e-07 3.677376e-07 0.9999998
[9,] 2.516904e-08 5.033808e-08 1.0000000
[10,] 3.359476e-09 6.718952e-09 1.0000000
[11,] 4.408008e-10 8.816015e-10 1.0000000
[12,] 5.285267e-11 1.057053e-10 1.0000000
[13,] 7.156520e-12 1.431304e-11 1.0000000
[14,] 8.215636e-13 1.643127e-12 1.0000000
[15,] 7.264401e-13 1.452880e-12 1.0000000
[16,] 3.550354e-13 7.100708e-13 1.0000000
[17,] 1.529644e-13 3.059288e-13 1.0000000
[18,] 5.124987e-13 1.024997e-12 1.0000000
[19,] 2.914435e-12 5.828871e-12 1.0000000
[20,] 3.540224e-10 7.080449e-10 1.0000000
[21,] 8.409323e-09 1.681865e-08 1.0000000
[22,] 1.120822e-07 2.241643e-07 0.9999999
[23,] 1.355845e-05 2.711691e-05 0.9999864
[24,] 1.157054e-03 2.314108e-03 0.9988429
[25,] 1.288694e-02 2.577387e-02 0.9871131
[26,] 3.465975e-01 6.931950e-01 0.6534025
[27,] 8.363881e-01 3.272238e-01 0.1636119
> postscript(file="/var/www/html/rcomp/tmp/1i1dy1261387531.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/29vpd1261387531.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/3ufqq1261387531.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/43o3d1261387531.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/5k8jc1261387531.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> qqnorm(mysum$resid, main='Residual Normal Q-Q Plot')
> qqline(mysum$resid)
> grid()
> dev.off()
null device
1
> (myerror <- as.ts(mysum$resid))
Time Series:
Start = 1
End = 60
Frequency = 1
1 2 3 4 5
-0.2509003167 -0.4690970574 -0.3599166129 -0.1672937982 0.5010011944
6 7 8 9 10
0.2093285423 0.3046724236 0.5744100236 0.7605411050 0.7380165421
11 12 13 14 15
0.5797550906 -0.0009985685 -0.2930400249 -0.0617280246 -0.3345149876
16 17 18 19 20
-0.4032029874 -0.9655634021 -0.5419892390 -0.2700877279 -0.4057599057
21 22 23 24 25
-0.4574972685 -0.3475631650 -0.2174649132 -0.0873990167 -0.6397687696
26 27 28 29 30
-0.7138665471 -0.5650143991 -0.4814078806 -0.5040965918 -0.4823256879
31 32 33 34 35
-1.0248502509 -0.9767517619 -0.7798011250 -0.9518344289 -0.8343589919
36 37 38 39 40
-0.2357692439 -0.4115813670 -0.0225650708 0.1794023367 0.4431730035
41 42 43 44 45
1.3952386628 1.3252051217 2.2574349292 2.8686474919 3.0782209437
46 47 48 49 50
2.7574996400 3.0633347803 2.9250386020 1.5952904782 1.2672567000
51 52 53 54 55
1.0800436630 0.6087316627 -0.4265798632 -0.5102187371 -1.2671693740
56 57 58 59 60
-2.0605458479 -2.6014636552 -2.1961185881 -2.5912659658 -2.6008717729
> postscript(file="/var/www/html/rcomp/tmp/66ncp1261387531.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> dum <- cbind(lag(myerror,k=1),myerror)
> dum
Time Series:
Start = 0
End = 60
Frequency = 1
lag(myerror, k = 1) myerror
0 -0.2509003167 NA
1 -0.4690970574 -0.2509003167
2 -0.3599166129 -0.4690970574
3 -0.1672937982 -0.3599166129
4 0.5010011944 -0.1672937982
5 0.2093285423 0.5010011944
6 0.3046724236 0.2093285423
7 0.5744100236 0.3046724236
8 0.7605411050 0.5744100236
9 0.7380165421 0.7605411050
10 0.5797550906 0.7380165421
11 -0.0009985685 0.5797550906
12 -0.2930400249 -0.0009985685
13 -0.0617280246 -0.2930400249
14 -0.3345149876 -0.0617280246
15 -0.4032029874 -0.3345149876
16 -0.9655634021 -0.4032029874
17 -0.5419892390 -0.9655634021
18 -0.2700877279 -0.5419892390
19 -0.4057599057 -0.2700877279
20 -0.4574972685 -0.4057599057
21 -0.3475631650 -0.4574972685
22 -0.2174649132 -0.3475631650
23 -0.0873990167 -0.2174649132
24 -0.6397687696 -0.0873990167
25 -0.7138665471 -0.6397687696
26 -0.5650143991 -0.7138665471
27 -0.4814078806 -0.5650143991
28 -0.5040965918 -0.4814078806
29 -0.4823256879 -0.5040965918
30 -1.0248502509 -0.4823256879
31 -0.9767517619 -1.0248502509
32 -0.7798011250 -0.9767517619
33 -0.9518344289 -0.7798011250
34 -0.8343589919 -0.9518344289
35 -0.2357692439 -0.8343589919
36 -0.4115813670 -0.2357692439
37 -0.0225650708 -0.4115813670
38 0.1794023367 -0.0225650708
39 0.4431730035 0.1794023367
40 1.3952386628 0.4431730035
41 1.3252051217 1.3952386628
42 2.2574349292 1.3252051217
43 2.8686474919 2.2574349292
44 3.0782209437 2.8686474919
45 2.7574996400 3.0782209437
46 3.0633347803 2.7574996400
47 2.9250386020 3.0633347803
48 1.5952904782 2.9250386020
49 1.2672567000 1.5952904782
50 1.0800436630 1.2672567000
51 0.6087316627 1.0800436630
52 -0.4265798632 0.6087316627
53 -0.5102187371 -0.4265798632
54 -1.2671693740 -0.5102187371
55 -2.0605458479 -1.2671693740
56 -2.6014636552 -2.0605458479
57 -2.1961185881 -2.6014636552
58 -2.5912659658 -2.1961185881
59 -2.6008717729 -2.5912659658
60 NA -2.6008717729
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -0.4690970574 -0.2509003167
[2,] -0.3599166129 -0.4690970574
[3,] -0.1672937982 -0.3599166129
[4,] 0.5010011944 -0.1672937982
[5,] 0.2093285423 0.5010011944
[6,] 0.3046724236 0.2093285423
[7,] 0.5744100236 0.3046724236
[8,] 0.7605411050 0.5744100236
[9,] 0.7380165421 0.7605411050
[10,] 0.5797550906 0.7380165421
[11,] -0.0009985685 0.5797550906
[12,] -0.2930400249 -0.0009985685
[13,] -0.0617280246 -0.2930400249
[14,] -0.3345149876 -0.0617280246
[15,] -0.4032029874 -0.3345149876
[16,] -0.9655634021 -0.4032029874
[17,] -0.5419892390 -0.9655634021
[18,] -0.2700877279 -0.5419892390
[19,] -0.4057599057 -0.2700877279
[20,] -0.4574972685 -0.4057599057
[21,] -0.3475631650 -0.4574972685
[22,] -0.2174649132 -0.3475631650
[23,] -0.0873990167 -0.2174649132
[24,] -0.6397687696 -0.0873990167
[25,] -0.7138665471 -0.6397687696
[26,] -0.5650143991 -0.7138665471
[27,] -0.4814078806 -0.5650143991
[28,] -0.5040965918 -0.4814078806
[29,] -0.4823256879 -0.5040965918
[30,] -1.0248502509 -0.4823256879
[31,] -0.9767517619 -1.0248502509
[32,] -0.7798011250 -0.9767517619
[33,] -0.9518344289 -0.7798011250
[34,] -0.8343589919 -0.9518344289
[35,] -0.2357692439 -0.8343589919
[36,] -0.4115813670 -0.2357692439
[37,] -0.0225650708 -0.4115813670
[38,] 0.1794023367 -0.0225650708
[39,] 0.4431730035 0.1794023367
[40,] 1.3952386628 0.4431730035
[41,] 1.3252051217 1.3952386628
[42,] 2.2574349292 1.3252051217
[43,] 2.8686474919 2.2574349292
[44,] 3.0782209437 2.8686474919
[45,] 2.7574996400 3.0782209437
[46,] 3.0633347803 2.7574996400
[47,] 2.9250386020 3.0633347803
[48,] 1.5952904782 2.9250386020
[49,] 1.2672567000 1.5952904782
[50,] 1.0800436630 1.2672567000
[51,] 0.6087316627 1.0800436630
[52,] -0.4265798632 0.6087316627
[53,] -0.5102187371 -0.4265798632
[54,] -1.2671693740 -0.5102187371
[55,] -2.0605458479 -1.2671693740
[56,] -2.6014636552 -2.0605458479
[57,] -2.1961185881 -2.6014636552
[58,] -2.5912659658 -2.1961185881
[59,] -2.6008717729 -2.5912659658
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -0.4690970574 -0.2509003167
2 -0.3599166129 -0.4690970574
3 -0.1672937982 -0.3599166129
4 0.5010011944 -0.1672937982
5 0.2093285423 0.5010011944
6 0.3046724236 0.2093285423
7 0.5744100236 0.3046724236
8 0.7605411050 0.5744100236
9 0.7380165421 0.7605411050
10 0.5797550906 0.7380165421
11 -0.0009985685 0.5797550906
12 -0.2930400249 -0.0009985685
13 -0.0617280246 -0.2930400249
14 -0.3345149876 -0.0617280246
15 -0.4032029874 -0.3345149876
16 -0.9655634021 -0.4032029874
17 -0.5419892390 -0.9655634021
18 -0.2700877279 -0.5419892390
19 -0.4057599057 -0.2700877279
20 -0.4574972685 -0.4057599057
21 -0.3475631650 -0.4574972685
22 -0.2174649132 -0.3475631650
23 -0.0873990167 -0.2174649132
24 -0.6397687696 -0.0873990167
25 -0.7138665471 -0.6397687696
26 -0.5650143991 -0.7138665471
27 -0.4814078806 -0.5650143991
28 -0.5040965918 -0.4814078806
29 -0.4823256879 -0.5040965918
30 -1.0248502509 -0.4823256879
31 -0.9767517619 -1.0248502509
32 -0.7798011250 -0.9767517619
33 -0.9518344289 -0.7798011250
34 -0.8343589919 -0.9518344289
35 -0.2357692439 -0.8343589919
36 -0.4115813670 -0.2357692439
37 -0.0225650708 -0.4115813670
38 0.1794023367 -0.0225650708
39 0.4431730035 0.1794023367
40 1.3952386628 0.4431730035
41 1.3252051217 1.3952386628
42 2.2574349292 1.3252051217
43 2.8686474919 2.2574349292
44 3.0782209437 2.8686474919
45 2.7574996400 3.0782209437
46 3.0633347803 2.7574996400
47 2.9250386020 3.0633347803
48 1.5952904782 2.9250386020
49 1.2672567000 1.5952904782
50 1.0800436630 1.2672567000
51 0.6087316627 1.0800436630
52 -0.4265798632 0.6087316627
53 -0.5102187371 -0.4265798632
54 -1.2671693740 -0.5102187371
55 -2.0605458479 -1.2671693740
56 -2.6014636552 -2.0605458479
57 -2.1961185881 -2.6014636552
58 -2.5912659658 -2.1961185881
59 -2.6008717729 -2.5912659658
> 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/7uxdv1261387532.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/8v6cx1261387532.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/9pebu1261387532.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/10pli01261387532.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/11oxxa1261387532.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/12ligy1261387532.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/13prq81261387532.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/1494um1261387532.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/154mtk1261387532.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/16dmah1261387532.tab")
+ }
>
> try(system("convert tmp/1i1dy1261387531.ps tmp/1i1dy1261387531.png",intern=TRUE))
character(0)
> try(system("convert tmp/29vpd1261387531.ps tmp/29vpd1261387531.png",intern=TRUE))
character(0)
> try(system("convert tmp/3ufqq1261387531.ps tmp/3ufqq1261387531.png",intern=TRUE))
character(0)
> try(system("convert tmp/43o3d1261387531.ps tmp/43o3d1261387531.png",intern=TRUE))
character(0)
> try(system("convert tmp/5k8jc1261387531.ps tmp/5k8jc1261387531.png",intern=TRUE))
character(0)
> try(system("convert tmp/66ncp1261387531.ps tmp/66ncp1261387531.png",intern=TRUE))
character(0)
> try(system("convert tmp/7uxdv1261387532.ps tmp/7uxdv1261387532.png",intern=TRUE))
character(0)
> try(system("convert tmp/8v6cx1261387532.ps tmp/8v6cx1261387532.png",intern=TRUE))
character(0)
> try(system("convert tmp/9pebu1261387532.ps tmp/9pebu1261387532.png",intern=TRUE))
character(0)
> try(system("convert tmp/10pli01261387532.ps tmp/10pli01261387532.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
2.380 1.543 3.625