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(269645,0,267037,0,258113,0,262813,0,267413,0,267366,0,264777,0,258863,0,254844,0,254868,0,277267,0,285351,0,286602,0,283042,0,276687,0,277915,0,277128,0,277103,0,275037,0,270150,0,267140,0,264993,0,287259,0,291186,0,292300,0,288186,0,281477,0,282656,0,280190,0,280408,0,276836,0,275216,0,274352,0,271311,0,289802,0,290726,0,292300,0,278506,0,269826,0,265861,0,269034,0,264176,0,255198,0,253353,0,246057,0,235372,0,258556,0,260993,0,254663,0,250643,0,243422,0,247105,0,248541,0,245039,0,237080,0,237085,0,225554,0,226839,0,247934,0,248333,1,246969,1,245098,1,246263,1,255765,1,264319,1,268347,1,273046,1,273963,1,267430,1,271993,1,292710,1,295881,1),dim=c(2,72),dimnames=list(c('Y','x'),1:72))
> y <- array(NA,dim=c(2,72),dimnames=list(c('Y','x'),1:72))
> 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 269645 0 1 0 0 0 0 0 0 0 0 0 0 1
2 267037 0 0 1 0 0 0 0 0 0 0 0 0 2
3 258113 0 0 0 1 0 0 0 0 0 0 0 0 3
4 262813 0 0 0 0 1 0 0 0 0 0 0 0 4
5 267413 0 0 0 0 0 1 0 0 0 0 0 0 5
6 267366 0 0 0 0 0 0 1 0 0 0 0 0 6
7 264777 0 0 0 0 0 0 0 1 0 0 0 0 7
8 258863 0 0 0 0 0 0 0 0 1 0 0 0 8
9 254844 0 0 0 0 0 0 0 0 0 1 0 0 9
10 254868 0 0 0 0 0 0 0 0 0 0 1 0 10
11 277267 0 0 0 0 0 0 0 0 0 0 0 1 11
12 285351 0 0 0 0 0 0 0 0 0 0 0 0 12
13 286602 0 1 0 0 0 0 0 0 0 0 0 0 13
14 283042 0 0 1 0 0 0 0 0 0 0 0 0 14
15 276687 0 0 0 1 0 0 0 0 0 0 0 0 15
16 277915 0 0 0 0 1 0 0 0 0 0 0 0 16
17 277128 0 0 0 0 0 1 0 0 0 0 0 0 17
18 277103 0 0 0 0 0 0 1 0 0 0 0 0 18
19 275037 0 0 0 0 0 0 0 1 0 0 0 0 19
20 270150 0 0 0 0 0 0 0 0 1 0 0 0 20
21 267140 0 0 0 0 0 0 0 0 0 1 0 0 21
22 264993 0 0 0 0 0 0 0 0 0 0 1 0 22
23 287259 0 0 0 0 0 0 0 0 0 0 0 1 23
24 291186 0 0 0 0 0 0 0 0 0 0 0 0 24
25 292300 0 1 0 0 0 0 0 0 0 0 0 0 25
26 288186 0 0 1 0 0 0 0 0 0 0 0 0 26
27 281477 0 0 0 1 0 0 0 0 0 0 0 0 27
28 282656 0 0 0 0 1 0 0 0 0 0 0 0 28
29 280190 0 0 0 0 0 1 0 0 0 0 0 0 29
30 280408 0 0 0 0 0 0 1 0 0 0 0 0 30
31 276836 0 0 0 0 0 0 0 1 0 0 0 0 31
32 275216 0 0 0 0 0 0 0 0 1 0 0 0 32
33 274352 0 0 0 0 0 0 0 0 0 1 0 0 33
34 271311 0 0 0 0 0 0 0 0 0 0 1 0 34
35 289802 0 0 0 0 0 0 0 0 0 0 0 1 35
36 290726 0 0 0 0 0 0 0 0 0 0 0 0 36
37 292300 0 1 0 0 0 0 0 0 0 0 0 0 37
38 278506 0 0 1 0 0 0 0 0 0 0 0 0 38
39 269826 0 0 0 1 0 0 0 0 0 0 0 0 39
40 265861 0 0 0 0 1 0 0 0 0 0 0 0 40
41 269034 0 0 0 0 0 1 0 0 0 0 0 0 41
42 264176 0 0 0 0 0 0 1 0 0 0 0 0 42
43 255198 0 0 0 0 0 0 0 1 0 0 0 0 43
44 253353 0 0 0 0 0 0 0 0 1 0 0 0 44
45 246057 0 0 0 0 0 0 0 0 0 1 0 0 45
46 235372 0 0 0 0 0 0 0 0 0 0 1 0 46
47 258556 0 0 0 0 0 0 0 0 0 0 0 1 47
48 260993 0 0 0 0 0 0 0 0 0 0 0 0 48
49 254663 0 1 0 0 0 0 0 0 0 0 0 0 49
50 250643 0 0 1 0 0 0 0 0 0 0 0 0 50
51 243422 0 0 0 1 0 0 0 0 0 0 0 0 51
52 247105 0 0 0 0 1 0 0 0 0 0 0 0 52
53 248541 0 0 0 0 0 1 0 0 0 0 0 0 53
54 245039 0 0 0 0 0 0 1 0 0 0 0 0 54
55 237080 0 0 0 0 0 0 0 1 0 0 0 0 55
56 237085 0 0 0 0 0 0 0 0 1 0 0 0 56
57 225554 0 0 0 0 0 0 0 0 0 1 0 0 57
58 226839 0 0 0 0 0 0 0 0 0 0 1 0 58
59 247934 0 0 0 0 0 0 0 0 0 0 0 1 59
60 248333 1 0 0 0 0 0 0 0 0 0 0 0 60
61 246969 1 1 0 0 0 0 0 0 0 0 0 0 61
62 245098 1 0 1 0 0 0 0 0 0 0 0 0 62
63 246263 1 0 0 1 0 0 0 0 0 0 0 0 63
64 255765 1 0 0 0 1 0 0 0 0 0 0 0 64
65 264319 1 0 0 0 0 1 0 0 0 0 0 0 65
66 268347 1 0 0 0 0 0 1 0 0 0 0 0 66
67 273046 1 0 0 0 0 0 0 1 0 0 0 0 67
68 273963 1 0 0 0 0 0 0 0 1 0 0 0 68
69 267430 1 0 0 0 0 0 0 0 0 1 0 0 69
70 271993 1 0 0 0 0 0 0 0 0 0 1 0 70
71 292710 1 0 0 0 0 0 0 0 0 0 0 1 71
72 295881 1 0 0 0 0 0 0 0 0 0 0 0 72
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) x M1 M2 M3 M4
293474.4 14540.2 -7702.2 -12230.6 -17885.2 -14697.9
M5 M6 M7 M8 M9 M10
-11813.5 -12045.1 -14989.8 -16747.7 -21823.8 -23024.5
M11 t
-1199.7 -466.1
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-31716 -11403 2534 12511 23773
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 293474.4 7110.3 41.274 < 2e-16 ***
x 14540.2 6151.3 2.364 0.021462 *
M1 -7702.2 8587.5 -0.897 0.373477
M2 -12230.6 8580.7 -1.425 0.159410
M3 -17885.2 8575.3 -2.086 0.041419 *
M4 -14697.9 8571.5 -1.715 0.091730 .
M5 -11813.5 8569.2 -1.379 0.173311
M6 -12045.1 8568.4 -1.406 0.165132
M7 -14989.8 8569.2 -1.749 0.085534 .
M8 -16747.7 8571.5 -1.954 0.055542 .
M9 -21823.8 8575.3 -2.545 0.013612 *
M10 -23024.5 8580.7 -2.683 0.009484 **
M11 -1199.7 8587.5 -0.140 0.889376
t -466.1 114.6 -4.066 0.000146 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 14780 on 58 degrees of freedom
Multiple R-squared: 0.3652, Adjusted R-squared: 0.2229
F-statistic: 2.567 on 13 and 58 DF, p-value: 0.007188
> 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.056439e-02 2.112878e-02 0.989435611
[2,] 4.129553e-03 8.259106e-03 0.995870447
[3,] 1.194970e-03 2.389940e-03 0.998805030
[4,] 2.825919e-04 5.651837e-04 0.999717408
[5,] 5.621489e-05 1.124298e-04 0.999943785
[6,] 1.541807e-05 3.083613e-05 0.999984582
[7,] 4.263773e-06 8.527545e-06 0.999995736
[8,] 3.669848e-06 7.339696e-06 0.999996330
[9,] 1.194377e-06 2.388753e-06 0.999998806
[10,] 4.264746e-07 8.529492e-07 0.999999574
[11,] 1.022109e-07 2.044217e-07 0.999999898
[12,] 3.377680e-08 6.755360e-08 0.999999966
[13,] 4.837717e-08 9.675434e-08 0.999999952
[14,] 3.458753e-08 6.917505e-08 0.999999965
[15,] 2.748608e-08 5.497215e-08 0.999999973
[16,] 7.124270e-09 1.424854e-08 0.999999993
[17,] 1.487133e-09 2.974267e-09 0.999999999
[18,] 3.291620e-10 6.583239e-10 1.000000000
[19,] 1.494806e-10 2.989612e-10 1.000000000
[20,] 4.721086e-10 9.442172e-10 1.000000000
[21,] 1.542280e-09 3.084560e-09 0.999999998
[22,] 1.175871e-07 2.351741e-07 0.999999882
[23,] 2.057087e-06 4.114173e-06 0.999997943
[24,] 3.904711e-05 7.809423e-05 0.999960953
[25,] 1.356159e-04 2.712318e-04 0.999864384
[26,] 5.685396e-04 1.137079e-03 0.999431460
[27,] 2.478906e-03 4.957812e-03 0.997521094
[28,] 4.596791e-03 9.193583e-03 0.995403209
[29,] 1.212302e-02 2.424605e-02 0.987876977
[30,] 3.010700e-02 6.021401e-02 0.969892996
[31,] 6.201653e-02 1.240331e-01 0.937983469
[32,] 1.637139e-01 3.274277e-01 0.836286137
[33,] 3.324465e-01 6.648930e-01 0.667553522
[34,] 5.710510e-01 8.578981e-01 0.428949049
[35,] 7.503224e-01 4.993552e-01 0.249677580
[36,] 8.831880e-01 2.336240e-01 0.116812022
[37,] 9.593975e-01 8.120494e-02 0.040602469
[38,] 9.939340e-01 1.213200e-02 0.006066001
[39,] 9.880399e-01 2.392026e-02 0.011960128
> postscript(file="/var/www/html/rcomp/tmp/10l351258740573.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/2ta501258740573.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/3ohwg1258740573.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/4176w1258740573.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/5s8bj1258740573.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 = 72
Frequency = 1
1 2 3 4 5 6 7
-15661.085 -13274.585 -16077.918 -14099.085 -11917.418 -11266.752 -10444.918
8 9 10 11 12 13 14
-14134.918 -12611.752 -10920.918 -9880.585 -2530.227 6889.092 8323.592
15 16 17 18 19 20 21
8089.259 6596.092 3390.759 4063.426 5408.259 2745.259 5277.426
22 23 24 25 26 27 28
4797.259 5704.592 8897.951 18180.270 19060.770 18472.437 16930.270
29 30 31 32 33 34 35
12045.937 12961.603 12800.437 13404.437 18082.603 16708.437 13840.770
36 37 38 39 40 41 42
14031.128 23773.447 14973.947 12414.614 5728.447 6483.114 2322.781
43 44 45 46 47 48 49
-3244.386 -2865.386 -4619.219 -13637.386 -11812.053 -10108.694 -8270.375
50 51 52 53 54 55 56
-7295.875 -8396.209 -7434.375 -8416.709 -11221.042 -15769.209 -13540.209
57 58 59 60 61 62 63
-19529.042 -16577.209 -16840.875 -31715.668 -24911.349 -21787.849 -14502.183
64 65 66 67 68 69 70
-7721.349 -1585.683 3139.984 11249.817 14390.817 13399.984 19629.817
71 72
18988.151 21425.509
> postscript(file="/var/www/html/rcomp/tmp/6ndkm1258740573.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 = 72
Frequency = 1
lag(myerror, k = 1) myerror
0 -15661.085 NA
1 -13274.585 -15661.085
2 -16077.918 -13274.585
3 -14099.085 -16077.918
4 -11917.418 -14099.085
5 -11266.752 -11917.418
6 -10444.918 -11266.752
7 -14134.918 -10444.918
8 -12611.752 -14134.918
9 -10920.918 -12611.752
10 -9880.585 -10920.918
11 -2530.227 -9880.585
12 6889.092 -2530.227
13 8323.592 6889.092
14 8089.259 8323.592
15 6596.092 8089.259
16 3390.759 6596.092
17 4063.426 3390.759
18 5408.259 4063.426
19 2745.259 5408.259
20 5277.426 2745.259
21 4797.259 5277.426
22 5704.592 4797.259
23 8897.951 5704.592
24 18180.270 8897.951
25 19060.770 18180.270
26 18472.437 19060.770
27 16930.270 18472.437
28 12045.937 16930.270
29 12961.603 12045.937
30 12800.437 12961.603
31 13404.437 12800.437
32 18082.603 13404.437
33 16708.437 18082.603
34 13840.770 16708.437
35 14031.128 13840.770
36 23773.447 14031.128
37 14973.947 23773.447
38 12414.614 14973.947
39 5728.447 12414.614
40 6483.114 5728.447
41 2322.781 6483.114
42 -3244.386 2322.781
43 -2865.386 -3244.386
44 -4619.219 -2865.386
45 -13637.386 -4619.219
46 -11812.053 -13637.386
47 -10108.694 -11812.053
48 -8270.375 -10108.694
49 -7295.875 -8270.375
50 -8396.209 -7295.875
51 -7434.375 -8396.209
52 -8416.709 -7434.375
53 -11221.042 -8416.709
54 -15769.209 -11221.042
55 -13540.209 -15769.209
56 -19529.042 -13540.209
57 -16577.209 -19529.042
58 -16840.875 -16577.209
59 -31715.668 -16840.875
60 -24911.349 -31715.668
61 -21787.849 -24911.349
62 -14502.183 -21787.849
63 -7721.349 -14502.183
64 -1585.683 -7721.349
65 3139.984 -1585.683
66 11249.817 3139.984
67 14390.817 11249.817
68 13399.984 14390.817
69 19629.817 13399.984
70 18988.151 19629.817
71 21425.509 18988.151
72 NA 21425.509
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -13274.585 -15661.085
[2,] -16077.918 -13274.585
[3,] -14099.085 -16077.918
[4,] -11917.418 -14099.085
[5,] -11266.752 -11917.418
[6,] -10444.918 -11266.752
[7,] -14134.918 -10444.918
[8,] -12611.752 -14134.918
[9,] -10920.918 -12611.752
[10,] -9880.585 -10920.918
[11,] -2530.227 -9880.585
[12,] 6889.092 -2530.227
[13,] 8323.592 6889.092
[14,] 8089.259 8323.592
[15,] 6596.092 8089.259
[16,] 3390.759 6596.092
[17,] 4063.426 3390.759
[18,] 5408.259 4063.426
[19,] 2745.259 5408.259
[20,] 5277.426 2745.259
[21,] 4797.259 5277.426
[22,] 5704.592 4797.259
[23,] 8897.951 5704.592
[24,] 18180.270 8897.951
[25,] 19060.770 18180.270
[26,] 18472.437 19060.770
[27,] 16930.270 18472.437
[28,] 12045.937 16930.270
[29,] 12961.603 12045.937
[30,] 12800.437 12961.603
[31,] 13404.437 12800.437
[32,] 18082.603 13404.437
[33,] 16708.437 18082.603
[34,] 13840.770 16708.437
[35,] 14031.128 13840.770
[36,] 23773.447 14031.128
[37,] 14973.947 23773.447
[38,] 12414.614 14973.947
[39,] 5728.447 12414.614
[40,] 6483.114 5728.447
[41,] 2322.781 6483.114
[42,] -3244.386 2322.781
[43,] -2865.386 -3244.386
[44,] -4619.219 -2865.386
[45,] -13637.386 -4619.219
[46,] -11812.053 -13637.386
[47,] -10108.694 -11812.053
[48,] -8270.375 -10108.694
[49,] -7295.875 -8270.375
[50,] -8396.209 -7295.875
[51,] -7434.375 -8396.209
[52,] -8416.709 -7434.375
[53,] -11221.042 -8416.709
[54,] -15769.209 -11221.042
[55,] -13540.209 -15769.209
[56,] -19529.042 -13540.209
[57,] -16577.209 -19529.042
[58,] -16840.875 -16577.209
[59,] -31715.668 -16840.875
[60,] -24911.349 -31715.668
[61,] -21787.849 -24911.349
[62,] -14502.183 -21787.849
[63,] -7721.349 -14502.183
[64,] -1585.683 -7721.349
[65,] 3139.984 -1585.683
[66,] 11249.817 3139.984
[67,] 14390.817 11249.817
[68,] 13399.984 14390.817
[69,] 19629.817 13399.984
[70,] 18988.151 19629.817
[71,] 21425.509 18988.151
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -13274.585 -15661.085
2 -16077.918 -13274.585
3 -14099.085 -16077.918
4 -11917.418 -14099.085
5 -11266.752 -11917.418
6 -10444.918 -11266.752
7 -14134.918 -10444.918
8 -12611.752 -14134.918
9 -10920.918 -12611.752
10 -9880.585 -10920.918
11 -2530.227 -9880.585
12 6889.092 -2530.227
13 8323.592 6889.092
14 8089.259 8323.592
15 6596.092 8089.259
16 3390.759 6596.092
17 4063.426 3390.759
18 5408.259 4063.426
19 2745.259 5408.259
20 5277.426 2745.259
21 4797.259 5277.426
22 5704.592 4797.259
23 8897.951 5704.592
24 18180.270 8897.951
25 19060.770 18180.270
26 18472.437 19060.770
27 16930.270 18472.437
28 12045.937 16930.270
29 12961.603 12045.937
30 12800.437 12961.603
31 13404.437 12800.437
32 18082.603 13404.437
33 16708.437 18082.603
34 13840.770 16708.437
35 14031.128 13840.770
36 23773.447 14031.128
37 14973.947 23773.447
38 12414.614 14973.947
39 5728.447 12414.614
40 6483.114 5728.447
41 2322.781 6483.114
42 -3244.386 2322.781
43 -2865.386 -3244.386
44 -4619.219 -2865.386
45 -13637.386 -4619.219
46 -11812.053 -13637.386
47 -10108.694 -11812.053
48 -8270.375 -10108.694
49 -7295.875 -8270.375
50 -8396.209 -7295.875
51 -7434.375 -8396.209
52 -8416.709 -7434.375
53 -11221.042 -8416.709
54 -15769.209 -11221.042
55 -13540.209 -15769.209
56 -19529.042 -13540.209
57 -16577.209 -19529.042
58 -16840.875 -16577.209
59 -31715.668 -16840.875
60 -24911.349 -31715.668
61 -21787.849 -24911.349
62 -14502.183 -21787.849
63 -7721.349 -14502.183
64 -1585.683 -7721.349
65 3139.984 -1585.683
66 11249.817 3139.984
67 14390.817 11249.817
68 13399.984 14390.817
69 19629.817 13399.984
70 18988.151 19629.817
71 21425.509 18988.151
> 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/75f4l1258740574.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/8lgmf1258740574.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/9dmdy1258740574.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/10osin1258740574.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/11aalg1258740574.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/12wwe71258740574.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/13fdf41258740574.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/146tja1258740574.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/15qunx1258740574.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/16amkj1258740574.tab")
+ }
>
> system("convert tmp/10l351258740573.ps tmp/10l351258740573.png")
> system("convert tmp/2ta501258740573.ps tmp/2ta501258740573.png")
> system("convert tmp/3ohwg1258740573.ps tmp/3ohwg1258740573.png")
> system("convert tmp/4176w1258740573.ps tmp/4176w1258740573.png")
> system("convert tmp/5s8bj1258740573.ps tmp/5s8bj1258740573.png")
> system("convert tmp/6ndkm1258740573.ps tmp/6ndkm1258740573.png")
> system("convert tmp/75f4l1258740574.ps tmp/75f4l1258740574.png")
> system("convert tmp/8lgmf1258740574.ps tmp/8lgmf1258740574.png")
> system("convert tmp/9dmdy1258740574.ps tmp/9dmdy1258740574.png")
> system("convert tmp/10osin1258740574.ps tmp/10osin1258740574.png")
>
>
> proc.time()
user system elapsed
2.578 1.601 2.995