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(269285,8.2,269829,8,270911,7.5,266844,6.8,271244,6.5,269907,6.6,271296,7.6,270157,8,271322,8.1,267179,7.7,264101,7.5,265518,7.6,269419,7.8,268714,7.8,272482,7.8,268351,7.5,268175,7.5,270674,7.1,272764,7.5,272599,7.5,270333,7.6,270846,7.7,270491,7.7,269160,7.9,274027,8.1,273784,8.2,276663,8.2,274525,8.2,271344,7.9,271115,7.3,270798,6.9,273911,6.6,273985,6.7,271917,6.9,273338,7,270601,7.1,273547,7.2,275363,7.1,281229,6.9,277793,7,279913,6.8,282500,6.4,280041,6.7,282166,6.6,290304,6.4,283519,6.3,287816,6.2,285226,6.5,287595,6.8,289741,6.8,289148,6.4,288301,6.1,290155,5.8,289648,6.1,288225,7.2,289351,7.3,294735,6.9,305333,6.1,309030,5.8,310215,6.2,321935,7.1,325734,7.7,320846,7.9,323023,7.7,319753,7.4,321753,7.5,320757,8,324479,8.1),dim=c(2,68),dimnames=list(c('Y','X'),1:68))
> y <- array(NA,dim=c(2,68),dimnames=list(c('Y','X'),1:68))
> 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 = '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
1 269285 8.2 1 0 0 0 0 0 0 0 0 0 0
2 269829 8.0 0 1 0 0 0 0 0 0 0 0 0
3 270911 7.5 0 0 1 0 0 0 0 0 0 0 0
4 266844 6.8 0 0 0 1 0 0 0 0 0 0 0
5 271244 6.5 0 0 0 0 1 0 0 0 0 0 0
6 269907 6.6 0 0 0 0 0 1 0 0 0 0 0
7 271296 7.6 0 0 0 0 0 0 1 0 0 0 0
8 270157 8.0 0 0 0 0 0 0 0 1 0 0 0
9 271322 8.1 0 0 0 0 0 0 0 0 1 0 0
10 267179 7.7 0 0 0 0 0 0 0 0 0 1 0
11 264101 7.5 0 0 0 0 0 0 0 0 0 0 1
12 265518 7.6 0 0 0 0 0 0 0 0 0 0 0
13 269419 7.8 1 0 0 0 0 0 0 0 0 0 0
14 268714 7.8 0 1 0 0 0 0 0 0 0 0 0
15 272482 7.8 0 0 1 0 0 0 0 0 0 0 0
16 268351 7.5 0 0 0 1 0 0 0 0 0 0 0
17 268175 7.5 0 0 0 0 1 0 0 0 0 0 0
18 270674 7.1 0 0 0 0 0 1 0 0 0 0 0
19 272764 7.5 0 0 0 0 0 0 1 0 0 0 0
20 272599 7.5 0 0 0 0 0 0 0 1 0 0 0
21 270333 7.6 0 0 0 0 0 0 0 0 1 0 0
22 270846 7.7 0 0 0 0 0 0 0 0 0 1 0
23 270491 7.7 0 0 0 0 0 0 0 0 0 0 1
24 269160 7.9 0 0 0 0 0 0 0 0 0 0 0
25 274027 8.1 1 0 0 0 0 0 0 0 0 0 0
26 273784 8.2 0 1 0 0 0 0 0 0 0 0 0
27 276663 8.2 0 0 1 0 0 0 0 0 0 0 0
28 274525 8.2 0 0 0 1 0 0 0 0 0 0 0
29 271344 7.9 0 0 0 0 1 0 0 0 0 0 0
30 271115 7.3 0 0 0 0 0 1 0 0 0 0 0
31 270798 6.9 0 0 0 0 0 0 1 0 0 0 0
32 273911 6.6 0 0 0 0 0 0 0 1 0 0 0
33 273985 6.7 0 0 0 0 0 0 0 0 1 0 0
34 271917 6.9 0 0 0 0 0 0 0 0 0 1 0
35 273338 7.0 0 0 0 0 0 0 0 0 0 0 1
36 270601 7.1 0 0 0 0 0 0 0 0 0 0 0
37 273547 7.2 1 0 0 0 0 0 0 0 0 0 0
38 275363 7.1 0 1 0 0 0 0 0 0 0 0 0
39 281229 6.9 0 0 1 0 0 0 0 0 0 0 0
40 277793 7.0 0 0 0 1 0 0 0 0 0 0 0
41 279913 6.8 0 0 0 0 1 0 0 0 0 0 0
42 282500 6.4 0 0 0 0 0 1 0 0 0 0 0
43 280041 6.7 0 0 0 0 0 0 1 0 0 0 0
44 282166 6.6 0 0 0 0 0 0 0 1 0 0 0
45 290304 6.4 0 0 0 0 0 0 0 0 1 0 0
46 283519 6.3 0 0 0 0 0 0 0 0 0 1 0
47 287816 6.2 0 0 0 0 0 0 0 0 0 0 1
48 285226 6.5 0 0 0 0 0 0 0 0 0 0 0
49 287595 6.8 1 0 0 0 0 0 0 0 0 0 0
50 289741 6.8 0 1 0 0 0 0 0 0 0 0 0
51 289148 6.4 0 0 1 0 0 0 0 0 0 0 0
52 288301 6.1 0 0 0 1 0 0 0 0 0 0 0
53 290155 5.8 0 0 0 0 1 0 0 0 0 0 0
54 289648 6.1 0 0 0 0 0 1 0 0 0 0 0
55 288225 7.2 0 0 0 0 0 0 1 0 0 0 0
56 289351 7.3 0 0 0 0 0 0 0 1 0 0 0
57 294735 6.9 0 0 0 0 0 0 0 0 1 0 0
58 305333 6.1 0 0 0 0 0 0 0 0 0 1 0
59 309030 5.8 0 0 0 0 0 0 0 0 0 0 1
60 310215 6.2 0 0 0 0 0 0 0 0 0 0 0
61 321935 7.1 1 0 0 0 0 0 0 0 0 0 0
62 325734 7.7 0 1 0 0 0 0 0 0 0 0 0
63 320846 7.9 0 0 1 0 0 0 0 0 0 0 0
64 323023 7.7 0 0 0 1 0 0 0 0 0 0 0
65 319753 7.4 0 0 0 0 1 0 0 0 0 0 0
66 321753 7.5 0 0 0 0 0 1 0 0 0 0 0
67 320757 8.0 0 0 0 0 0 0 1 0 0 0 0
68 324479 8.1 0 0 0 0 0 0 0 1 0 0 0
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X M1 M2 M3 M4
317487.9 -5289.5 4994.4 6573.2 7132.1 3824.2
M5 M6 M7 M8 M9 M10
2881.1 2923.2 5193.8 6833.8 415.0 -1019.9
M11
-352.5
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-18499 -11153 -6656 1766 43002
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 317487.9 28089.2 11.303 5.79e-16 ***
X -5289.5 3793.8 -1.394 0.169
M1 4994.4 11596.3 0.431 0.668
M2 6573.2 11638.2 0.565 0.575
M3 7132.1 11551.6 0.617 0.540
M4 3824.2 11471.9 0.333 0.740
M5 2881.1 11460.1 0.251 0.802
M6 2923.2 11488.7 0.254 0.800
M7 5193.8 11497.8 0.452 0.653
M8 6833.8 11509.2 0.594 0.555
M9 415.0 11969.7 0.035 0.972
M10 -1019.9 11974.5 -0.085 0.932
M11 -352.5 11994.9 -0.029 0.977
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 18920 on 55 degrees of freedom
Multiple R-squared: 0.04545, Adjusted R-squared: -0.1628
F-statistic: 0.2182 on 12 and 55 DF, p-value: 0.9969
> if (n > n25) {
+ kp3 <- k + 3
+ nmkm3 <- n - k - 3
+ gqarr <- array(NA, dim=c(nmkm3-kp3+1,3))
+ numgqtests <- 0
+ numsignificant1 <- 0
+ numsignificant5 <- 0
+ numsignificant10 <- 0
+ for (mypoint in kp3:nmkm3) {
+ j <- 0
+ numgqtests <- numgqtests + 1
+ for (myalt in c('greater', 'two.sided', 'less')) {
+ j <- j + 1
+ gqarr[mypoint-kp3+1,j] <- gqtest(mylm, point=mypoint, alternative=myalt)$p.value
+ }
+ if (gqarr[mypoint-kp3+1,2] < 0.01) numsignificant1 <- numsignificant1 + 1
+ if (gqarr[mypoint-kp3+1,2] < 0.05) numsignificant5 <- numsignificant5 + 1
+ if (gqarr[mypoint-kp3+1,2] < 0.10) numsignificant10 <- numsignificant10 + 1
+ }
+ gqarr
+ }
[,1] [,2] [,3]
[1,] 4.047264e-05 8.094528e-05 0.9999595
[2,] 3.963516e-05 7.927033e-05 0.9999604
[3,] 2.789908e-06 5.579817e-06 0.9999972
[4,] 2.199791e-07 4.399581e-07 0.9999998
[5,] 2.562896e-08 5.125791e-08 1.0000000
[6,] 2.062490e-09 4.124980e-09 1.0000000
[7,] 8.604355e-10 1.720871e-09 1.0000000
[8,] 2.564418e-09 5.128836e-09 1.0000000
[9,] 5.967867e-10 1.193573e-09 1.0000000
[10,] 3.210632e-10 6.421265e-10 1.0000000
[11,] 1.379508e-10 2.759016e-10 1.0000000
[12,] 5.376750e-11 1.075350e-10 1.0000000
[13,] 2.352996e-11 4.705992e-11 1.0000000
[14,] 7.354483e-12 1.470897e-11 1.0000000
[15,] 2.070984e-12 4.141969e-12 1.0000000
[16,] 3.276823e-13 6.553646e-13 1.0000000
[17,] 2.249426e-13 4.498853e-13 1.0000000
[18,] 1.229133e-13 2.458267e-13 1.0000000
[19,] 7.511781e-14 1.502356e-13 1.0000000
[20,] 4.030467e-13 8.060935e-13 1.0000000
[21,] 9.987605e-13 1.997521e-12 1.0000000
[22,] 1.938464e-12 3.876927e-12 1.0000000
[23,] 3.786484e-12 7.572968e-12 1.0000000
[24,] 1.080725e-11 2.161450e-11 1.0000000
[25,] 1.357262e-10 2.714523e-10 1.0000000
[26,] 4.674348e-09 9.348697e-09 1.0000000
[27,] 3.315220e-08 6.630441e-08 1.0000000
[28,] 2.340342e-08 4.680684e-08 1.0000000
[29,] 1.816225e-08 3.632450e-08 1.0000000
[30,] 1.311418e-07 2.622835e-07 0.9999999
[31,] 6.708958e-07 1.341792e-06 0.9999993
[32,] 2.367434e-05 4.734868e-05 0.9999763
[33,] 1.651794e-03 3.303588e-03 0.9983482
[34,] 4.430116e-02 8.860233e-02 0.9556988
[35,] 1.022510e-01 2.045020e-01 0.8977490
[36,] 6.259830e-02 1.251966e-01 0.9374017
[37,] 3.648449e-02 7.296898e-02 0.9635155
> postscript(file="/var/www/html/rcomp/tmp/1isa31258805121.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/2gmu11258805121.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/3bxoy1258805121.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/4fw2z1258805121.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/5ar0r1258805121.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 = 68
Frequency = 1
1 2 3 4 5 6
-9823.3300 -11916.0313 -14037.6914 -18499.4604 -14743.2608 -15593.3845
7 8 9 10 11 12
-11185.4736 -11848.6551 -3735.8752 -8559.7762 -13363.1267 -11769.6673
13 14 15 16 17 18
-11805.1320 -14088.9323 -10879.8399 -13289.8069 -12522.7557 -12181.6320
19 20 21 22 23 24
-10246.4241 -12051.4076 -7369.6277 -4892.7762 -5915.2257 -6540.8158
25 26 27 28 29 30
-5610.2805 -6903.1303 -4583.0379 -3413.1534 -7237.9537 -10682.7310
31 32 33 34 35 36
-15386.1271 -15499.9621 -8478.1822 -8053.3802 -6770.8792 -9331.4198
37 38 39 40 41 42
-10850.8350 -11142.5858 -6893.3944 -6492.5594 -4487.4093 -4058.2855
43 44 45 46 47 48
-7201.0281 -7244.9621 6253.9663 374.9168 3475.5168 2119.8772
49 50 51 52 53 54
1081.3630 1648.5627 -1619.1469 -745.1139 465.0857 1502.8630
55 56 57 58 59 60
3627.7244 3642.6914 13329.7188 21131.0158 22573.7148 25522.0257
61 62 63 64 65 66
37008.2145 42402.1172 38013.1106 42440.0941 38526.2938 41013.1700
67 68
40391.3284 43002.2954
> postscript(file="/var/www/html/rcomp/tmp/60v7u1258805121.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 = 68
Frequency = 1
lag(myerror, k = 1) myerror
0 -9823.3300 NA
1 -11916.0313 -9823.3300
2 -14037.6914 -11916.0313
3 -18499.4604 -14037.6914
4 -14743.2608 -18499.4604
5 -15593.3845 -14743.2608
6 -11185.4736 -15593.3845
7 -11848.6551 -11185.4736
8 -3735.8752 -11848.6551
9 -8559.7762 -3735.8752
10 -13363.1267 -8559.7762
11 -11769.6673 -13363.1267
12 -11805.1320 -11769.6673
13 -14088.9323 -11805.1320
14 -10879.8399 -14088.9323
15 -13289.8069 -10879.8399
16 -12522.7557 -13289.8069
17 -12181.6320 -12522.7557
18 -10246.4241 -12181.6320
19 -12051.4076 -10246.4241
20 -7369.6277 -12051.4076
21 -4892.7762 -7369.6277
22 -5915.2257 -4892.7762
23 -6540.8158 -5915.2257
24 -5610.2805 -6540.8158
25 -6903.1303 -5610.2805
26 -4583.0379 -6903.1303
27 -3413.1534 -4583.0379
28 -7237.9537 -3413.1534
29 -10682.7310 -7237.9537
30 -15386.1271 -10682.7310
31 -15499.9621 -15386.1271
32 -8478.1822 -15499.9621
33 -8053.3802 -8478.1822
34 -6770.8792 -8053.3802
35 -9331.4198 -6770.8792
36 -10850.8350 -9331.4198
37 -11142.5858 -10850.8350
38 -6893.3944 -11142.5858
39 -6492.5594 -6893.3944
40 -4487.4093 -6492.5594
41 -4058.2855 -4487.4093
42 -7201.0281 -4058.2855
43 -7244.9621 -7201.0281
44 6253.9663 -7244.9621
45 374.9168 6253.9663
46 3475.5168 374.9168
47 2119.8772 3475.5168
48 1081.3630 2119.8772
49 1648.5627 1081.3630
50 -1619.1469 1648.5627
51 -745.1139 -1619.1469
52 465.0857 -745.1139
53 1502.8630 465.0857
54 3627.7244 1502.8630
55 3642.6914 3627.7244
56 13329.7188 3642.6914
57 21131.0158 13329.7188
58 22573.7148 21131.0158
59 25522.0257 22573.7148
60 37008.2145 25522.0257
61 42402.1172 37008.2145
62 38013.1106 42402.1172
63 42440.0941 38013.1106
64 38526.2938 42440.0941
65 41013.1700 38526.2938
66 40391.3284 41013.1700
67 43002.2954 40391.3284
68 NA 43002.2954
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -11916.0313 -9823.3300
[2,] -14037.6914 -11916.0313
[3,] -18499.4604 -14037.6914
[4,] -14743.2608 -18499.4604
[5,] -15593.3845 -14743.2608
[6,] -11185.4736 -15593.3845
[7,] -11848.6551 -11185.4736
[8,] -3735.8752 -11848.6551
[9,] -8559.7762 -3735.8752
[10,] -13363.1267 -8559.7762
[11,] -11769.6673 -13363.1267
[12,] -11805.1320 -11769.6673
[13,] -14088.9323 -11805.1320
[14,] -10879.8399 -14088.9323
[15,] -13289.8069 -10879.8399
[16,] -12522.7557 -13289.8069
[17,] -12181.6320 -12522.7557
[18,] -10246.4241 -12181.6320
[19,] -12051.4076 -10246.4241
[20,] -7369.6277 -12051.4076
[21,] -4892.7762 -7369.6277
[22,] -5915.2257 -4892.7762
[23,] -6540.8158 -5915.2257
[24,] -5610.2805 -6540.8158
[25,] -6903.1303 -5610.2805
[26,] -4583.0379 -6903.1303
[27,] -3413.1534 -4583.0379
[28,] -7237.9537 -3413.1534
[29,] -10682.7310 -7237.9537
[30,] -15386.1271 -10682.7310
[31,] -15499.9621 -15386.1271
[32,] -8478.1822 -15499.9621
[33,] -8053.3802 -8478.1822
[34,] -6770.8792 -8053.3802
[35,] -9331.4198 -6770.8792
[36,] -10850.8350 -9331.4198
[37,] -11142.5858 -10850.8350
[38,] -6893.3944 -11142.5858
[39,] -6492.5594 -6893.3944
[40,] -4487.4093 -6492.5594
[41,] -4058.2855 -4487.4093
[42,] -7201.0281 -4058.2855
[43,] -7244.9621 -7201.0281
[44,] 6253.9663 -7244.9621
[45,] 374.9168 6253.9663
[46,] 3475.5168 374.9168
[47,] 2119.8772 3475.5168
[48,] 1081.3630 2119.8772
[49,] 1648.5627 1081.3630
[50,] -1619.1469 1648.5627
[51,] -745.1139 -1619.1469
[52,] 465.0857 -745.1139
[53,] 1502.8630 465.0857
[54,] 3627.7244 1502.8630
[55,] 3642.6914 3627.7244
[56,] 13329.7188 3642.6914
[57,] 21131.0158 13329.7188
[58,] 22573.7148 21131.0158
[59,] 25522.0257 22573.7148
[60,] 37008.2145 25522.0257
[61,] 42402.1172 37008.2145
[62,] 38013.1106 42402.1172
[63,] 42440.0941 38013.1106
[64,] 38526.2938 42440.0941
[65,] 41013.1700 38526.2938
[66,] 40391.3284 41013.1700
[67,] 43002.2954 40391.3284
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -11916.0313 -9823.3300
2 -14037.6914 -11916.0313
3 -18499.4604 -14037.6914
4 -14743.2608 -18499.4604
5 -15593.3845 -14743.2608
6 -11185.4736 -15593.3845
7 -11848.6551 -11185.4736
8 -3735.8752 -11848.6551
9 -8559.7762 -3735.8752
10 -13363.1267 -8559.7762
11 -11769.6673 -13363.1267
12 -11805.1320 -11769.6673
13 -14088.9323 -11805.1320
14 -10879.8399 -14088.9323
15 -13289.8069 -10879.8399
16 -12522.7557 -13289.8069
17 -12181.6320 -12522.7557
18 -10246.4241 -12181.6320
19 -12051.4076 -10246.4241
20 -7369.6277 -12051.4076
21 -4892.7762 -7369.6277
22 -5915.2257 -4892.7762
23 -6540.8158 -5915.2257
24 -5610.2805 -6540.8158
25 -6903.1303 -5610.2805
26 -4583.0379 -6903.1303
27 -3413.1534 -4583.0379
28 -7237.9537 -3413.1534
29 -10682.7310 -7237.9537
30 -15386.1271 -10682.7310
31 -15499.9621 -15386.1271
32 -8478.1822 -15499.9621
33 -8053.3802 -8478.1822
34 -6770.8792 -8053.3802
35 -9331.4198 -6770.8792
36 -10850.8350 -9331.4198
37 -11142.5858 -10850.8350
38 -6893.3944 -11142.5858
39 -6492.5594 -6893.3944
40 -4487.4093 -6492.5594
41 -4058.2855 -4487.4093
42 -7201.0281 -4058.2855
43 -7244.9621 -7201.0281
44 6253.9663 -7244.9621
45 374.9168 6253.9663
46 3475.5168 374.9168
47 2119.8772 3475.5168
48 1081.3630 2119.8772
49 1648.5627 1081.3630
50 -1619.1469 1648.5627
51 -745.1139 -1619.1469
52 465.0857 -745.1139
53 1502.8630 465.0857
54 3627.7244 1502.8630
55 3642.6914 3627.7244
56 13329.7188 3642.6914
57 21131.0158 13329.7188
58 22573.7148 21131.0158
59 25522.0257 22573.7148
60 37008.2145 25522.0257
61 42402.1172 37008.2145
62 38013.1106 42402.1172
63 42440.0941 38013.1106
64 38526.2938 42440.0941
65 41013.1700 38526.2938
66 40391.3284 41013.1700
67 43002.2954 40391.3284
> 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/761i21258805121.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/8nph51258805121.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/94fza1258805121.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/10elfx1258805121.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/11r6ym1258805121.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/128azm1258805121.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/1329xy1258805121.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/1464xx1258805121.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/158x8u1258805121.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/16hydm1258805121.tab")
+ }
>
> system("convert tmp/1isa31258805121.ps tmp/1isa31258805121.png")
> system("convert tmp/2gmu11258805121.ps tmp/2gmu11258805121.png")
> system("convert tmp/3bxoy1258805121.ps tmp/3bxoy1258805121.png")
> system("convert tmp/4fw2z1258805121.ps tmp/4fw2z1258805121.png")
> system("convert tmp/5ar0r1258805121.ps tmp/5ar0r1258805121.png")
> system("convert tmp/60v7u1258805121.ps tmp/60v7u1258805121.png")
> system("convert tmp/761i21258805121.ps tmp/761i21258805121.png")
> system("convert tmp/8nph51258805121.ps tmp/8nph51258805121.png")
> system("convert tmp/94fza1258805121.ps tmp/94fza1258805121.png")
> system("convert tmp/10elfx1258805121.ps tmp/10elfx1258805121.png")
>
>
> proc.time()
user system elapsed
2.516 1.571 2.919