R version 2.8.0 (2008-10-20)
Copyright (C) 2008 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(78.4,0,114.6,0,113.3,0,117.0,0,99.6,0,99.4,0,101.9,0,115.2,0,108.5,0,113.8,0,121.0,0,92.2,0,90.2,0,101.5,0,126.6,0,93.9,0,89.8,0,93.4,0,101.5,0,110.4,0,105.9,0,108.4,0,113.9,0,86.1,0,69.4,0,101.2,0,100.5,0,98.0,0,106.6,0,90.1,0,96.9,0,125.9,0,112.0,0,100.0,0,123.9,0,79.8,0,83.4,0,113.6,0,112.9,0,104.0,0,109.9,0,99.0,0,106.3,0,128.9,0,111.1,0,102.9,0,130.0,0,87.0,0,87.5,0,117.6,0,103.4,0,110.8,0,112.6,0,102.5,0,112.4,0,135.6,0,105.1,0,127.7,0,137.0,0,91.0,0,90.5,0,122.4,0,123.3,0,124.3,0,120.0,0,118.1,0,119.0,0,142.7,0,123.6,0,129.6,0,151.6,0,110.4,1,99.2,1,130.5,1,136.2,1,129.7,1,128.0,1,121.6,1,135.8,1,143.8,1,147.5,1,136.2,1,156.6,1,123.3,1,100.4,1),dim=c(2,85),dimnames=list(c('I','D'),1:85))
> y <- array(NA,dim=c(2,85),dimnames=list(c('I','D'),1:85))
> 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
I D M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 78.4 0 1 0 0 0 0 0 0 0 0 0 0 1
2 114.6 0 0 1 0 0 0 0 0 0 0 0 0 2
3 113.3 0 0 0 1 0 0 0 0 0 0 0 0 3
4 117.0 0 0 0 0 1 0 0 0 0 0 0 0 4
5 99.6 0 0 0 0 0 1 0 0 0 0 0 0 5
6 99.4 0 0 0 0 0 0 1 0 0 0 0 0 6
7 101.9 0 0 0 0 0 0 0 1 0 0 0 0 7
8 115.2 0 0 0 0 0 0 0 0 1 0 0 0 8
9 108.5 0 0 0 0 0 0 0 0 0 1 0 0 9
10 113.8 0 0 0 0 0 0 0 0 0 0 1 0 10
11 121.0 0 0 0 0 0 0 0 0 0 0 0 1 11
12 92.2 0 0 0 0 0 0 0 0 0 0 0 0 12
13 90.2 0 1 0 0 0 0 0 0 0 0 0 0 13
14 101.5 0 0 1 0 0 0 0 0 0 0 0 0 14
15 126.6 0 0 0 1 0 0 0 0 0 0 0 0 15
16 93.9 0 0 0 0 1 0 0 0 0 0 0 0 16
17 89.8 0 0 0 0 0 1 0 0 0 0 0 0 17
18 93.4 0 0 0 0 0 0 1 0 0 0 0 0 18
19 101.5 0 0 0 0 0 0 0 1 0 0 0 0 19
20 110.4 0 0 0 0 0 0 0 0 1 0 0 0 20
21 105.9 0 0 0 0 0 0 0 0 0 1 0 0 21
22 108.4 0 0 0 0 0 0 0 0 0 0 1 0 22
23 113.9 0 0 0 0 0 0 0 0 0 0 0 1 23
24 86.1 0 0 0 0 0 0 0 0 0 0 0 0 24
25 69.4 0 1 0 0 0 0 0 0 0 0 0 0 25
26 101.2 0 0 1 0 0 0 0 0 0 0 0 0 26
27 100.5 0 0 0 1 0 0 0 0 0 0 0 0 27
28 98.0 0 0 0 0 1 0 0 0 0 0 0 0 28
29 106.6 0 0 0 0 0 1 0 0 0 0 0 0 29
30 90.1 0 0 0 0 0 0 1 0 0 0 0 0 30
31 96.9 0 0 0 0 0 0 0 1 0 0 0 0 31
32 125.9 0 0 0 0 0 0 0 0 1 0 0 0 32
33 112.0 0 0 0 0 0 0 0 0 0 1 0 0 33
34 100.0 0 0 0 0 0 0 0 0 0 0 1 0 34
35 123.9 0 0 0 0 0 0 0 0 0 0 0 1 35
36 79.8 0 0 0 0 0 0 0 0 0 0 0 0 36
37 83.4 0 1 0 0 0 0 0 0 0 0 0 0 37
38 113.6 0 0 1 0 0 0 0 0 0 0 0 0 38
39 112.9 0 0 0 1 0 0 0 0 0 0 0 0 39
40 104.0 0 0 0 0 1 0 0 0 0 0 0 0 40
41 109.9 0 0 0 0 0 1 0 0 0 0 0 0 41
42 99.0 0 0 0 0 0 0 1 0 0 0 0 0 42
43 106.3 0 0 0 0 0 0 0 1 0 0 0 0 43
44 128.9 0 0 0 0 0 0 0 0 1 0 0 0 44
45 111.1 0 0 0 0 0 0 0 0 0 1 0 0 45
46 102.9 0 0 0 0 0 0 0 0 0 0 1 0 46
47 130.0 0 0 0 0 0 0 0 0 0 0 0 1 47
48 87.0 0 0 0 0 0 0 0 0 0 0 0 0 48
49 87.5 0 1 0 0 0 0 0 0 0 0 0 0 49
50 117.6 0 0 1 0 0 0 0 0 0 0 0 0 50
51 103.4 0 0 0 1 0 0 0 0 0 0 0 0 51
52 110.8 0 0 0 0 1 0 0 0 0 0 0 0 52
53 112.6 0 0 0 0 0 1 0 0 0 0 0 0 53
54 102.5 0 0 0 0 0 0 1 0 0 0 0 0 54
55 112.4 0 0 0 0 0 0 0 1 0 0 0 0 55
56 135.6 0 0 0 0 0 0 0 0 1 0 0 0 56
57 105.1 0 0 0 0 0 0 0 0 0 1 0 0 57
58 127.7 0 0 0 0 0 0 0 0 0 0 1 0 58
59 137.0 0 0 0 0 0 0 0 0 0 0 0 1 59
60 91.0 0 0 0 0 0 0 0 0 0 0 0 0 60
61 90.5 0 1 0 0 0 0 0 0 0 0 0 0 61
62 122.4 0 0 1 0 0 0 0 0 0 0 0 0 62
63 123.3 0 0 0 1 0 0 0 0 0 0 0 0 63
64 124.3 0 0 0 0 1 0 0 0 0 0 0 0 64
65 120.0 0 0 0 0 0 1 0 0 0 0 0 0 65
66 118.1 0 0 0 0 0 0 1 0 0 0 0 0 66
67 119.0 0 0 0 0 0 0 0 1 0 0 0 0 67
68 142.7 0 0 0 0 0 0 0 0 1 0 0 0 68
69 123.6 0 0 0 0 0 0 0 0 0 1 0 0 69
70 129.6 0 0 0 0 0 0 0 0 0 0 1 0 70
71 151.6 0 0 0 0 0 0 0 0 0 0 0 1 71
72 110.4 1 0 0 0 0 0 0 0 0 0 0 0 72
73 99.2 1 1 0 0 0 0 0 0 0 0 0 0 73
74 130.5 1 0 1 0 0 0 0 0 0 0 0 0 74
75 136.2 1 0 0 1 0 0 0 0 0 0 0 0 75
76 129.7 1 0 0 0 1 0 0 0 0 0 0 0 76
77 128.0 1 0 0 0 0 1 0 0 0 0 0 0 77
78 121.6 1 0 0 0 0 0 1 0 0 0 0 0 78
79 135.8 1 0 0 0 0 0 0 1 0 0 0 0 79
80 143.8 1 0 0 0 0 0 0 0 1 0 0 0 80
81 147.5 1 0 0 0 0 0 0 0 0 1 0 0 81
82 136.2 1 0 0 0 0 0 0 0 0 0 1 0 82
83 156.6 1 0 0 0 0 0 0 0 0 0 0 1 83
84 123.3 1 0 0 0 0 0 0 0 0 0 0 0 84
85 100.4 1 1 0 0 0 0 0 0 0 0 0 0 85
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) D M1 M2 M3 M4
79.1306 11.9203 -6.5153 23.2423 25.0827 19.3087
M5 M6 M7 M8 M9 M10
17.4348 11.1037 17.9298 36.0415 23.0819 23.5079
M11 t
39.7197 0.2739
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-14.7844 -4.5672 -0.2653 3.5854 18.2776
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 79.13057 3.41207 23.191 < 2e-16 ***
D 11.92031 2.99682 3.978 0.000166 ***
M1 -6.51527 4.01772 -1.622 0.109315
M2 23.24235 4.16203 5.584 4.04e-07 ***
M3 25.08269 4.16036 6.029 6.69e-08 ***
M4 19.30874 4.15918 4.642 1.54e-05 ***
M5 17.43480 4.15850 4.193 7.84e-05 ***
M6 11.10371 4.15831 2.670 0.009390 **
M7 17.92977 4.15861 4.311 5.14e-05 ***
M8 36.04154 4.15940 8.665 9.78e-13 ***
M9 23.08188 4.16069 5.548 4.68e-07 ***
M10 23.50793 4.16247 5.648 3.14e-07 ***
M11 39.71970 4.16474 9.537 2.41e-14 ***
t 0.27394 0.04527 6.052 6.10e-08 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 7.755 on 71 degrees of freedom
Multiple R-squared: 0.8351, Adjusted R-squared: 0.8049
F-statistic: 27.65 on 13 and 71 DF, p-value: < 2.2e-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.9978907 0.004218601 0.002109301
[2,] 0.9950339 0.009932114 0.004966057
[3,] 0.9908159 0.018368166 0.009184083
[4,] 0.9820707 0.035858516 0.017929258
[5,] 0.9695268 0.060946354 0.030473177
[6,] 0.9587663 0.082467422 0.041233711
[7,] 0.9377101 0.124579711 0.062289855
[8,] 0.9236021 0.152795756 0.076397878
[9,] 0.9109775 0.178045033 0.089022516
[10,] 0.8702062 0.259587507 0.129793753
[11,] 0.8838740 0.232251957 0.116125979
[12,] 0.8389329 0.322134228 0.161067114
[13,] 0.9573263 0.085347460 0.042673730
[14,] 0.9357737 0.128452540 0.064226270
[15,] 0.9100143 0.179971322 0.089985661
[16,] 0.9651437 0.069712554 0.034856277
[17,] 0.9712365 0.057526911 0.028763456
[18,] 0.9643545 0.071290989 0.035645494
[19,] 0.9618804 0.076239109 0.038119555
[20,] 0.9489636 0.102072762 0.051036381
[21,] 0.9575206 0.084958899 0.042479449
[22,] 0.9642710 0.071458036 0.035729018
[23,] 0.9619978 0.076004393 0.038002196
[24,] 0.9458052 0.108389646 0.054194823
[25,] 0.9597372 0.080525543 0.040262771
[26,] 0.9475160 0.104968011 0.052484005
[27,] 0.9332828 0.133434354 0.066717177
[28,] 0.9423984 0.115203289 0.057601645
[29,] 0.9269633 0.146073432 0.073036716
[30,] 0.9259095 0.148180919 0.074090459
[31,] 0.9122653 0.175469407 0.087734703
[32,] 0.8789465 0.242106954 0.121053477
[33,] 0.8970672 0.205865537 0.102932768
[34,] 0.8914341 0.217131765 0.108565882
[35,] 0.9229259 0.154148167 0.077074084
[36,] 0.8958199 0.208360276 0.104180138
[37,] 0.8754507 0.249098689 0.124549344
[38,] 0.8346368 0.330726304 0.165363152
[39,] 0.7911043 0.417791346 0.208895673
[40,] 0.8030920 0.393816032 0.196908016
[41,] 0.9158749 0.168250132 0.084125066
[42,] 0.9423722 0.115255643 0.057627822
[43,] 0.9203707 0.159258630 0.079629315
[44,] 0.9606237 0.078752548 0.039376274
[45,] 0.9365744 0.126851115 0.063425557
[46,] 0.8978336 0.204332817 0.102166408
[47,] 0.8502107 0.299578632 0.149789316
[48,] 0.7995058 0.400988468 0.200494234
[49,] 0.7086565 0.582687027 0.291343514
[50,] 0.6537651 0.692469887 0.346234943
[51,] 0.5688773 0.862245448 0.431122724
[52,] 0.5268845 0.946230965 0.473115482
> postscript(file="/var/www/html/rcomp/tmp/1w2a31227728112.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/2gm5g1227728112.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/3n3ho1227728112.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/4hymi1227728112.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/5olhb1227728112.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 = 85
Frequency = 1
1 2 3 4 5 6
5.5107509 11.6791931 8.2649074 17.4649074 1.6649074 7.5220502
7 8 9 10 11 12
2.9220502 -2.1636641 3.8220502 8.4220502 -0.8636641 9.7820941
13 14 15 16 17 18
14.0234154 -4.7081424 18.2775719 -8.9224281 -11.4224281 -1.7652852
19 20 21 22 23 24
-0.7652852 -10.2509995 -2.0652852 -0.2652852 -11.2509995 0.3947587
25 26 27 28 29 30
-10.0639200 -8.2954778 -11.1097635 -8.1097635 2.0902365 -8.3526207
31 32 33 34 35 36
-8.6526207 1.9616650 0.7473793 -11.9526207 -4.5383350 -9.1925768
37 38 39 40 41 42
0.6487445 0.8171867 -1.9970990 -5.3970990 2.1029010 -2.7399561
43 44 45 46 47 48
-2.5399561 1.6743296 -3.4399561 -12.3399561 -1.7256704 -5.2799122
49 50 51 52 53 54
1.4614091 1.5298513 -14.7844344 -1.8844344 1.5155656 -2.5272916
55 56 57 58 59 60
0.2727084 5.0869941 -12.7272916 9.1727084 1.9869941 -4.5672477
61 62 63 64 65 66
1.1740736 3.0425158 1.8282301 8.3282301 5.6282301 9.7853730
67 68 69 70 71 72
3.5853730 8.8996587 2.4853730 7.7853730 13.2996587 -0.3748903
73 74 75 76 77 78
-5.3335690 -4.0651268 -0.4794125 -1.4794125 -1.5794125 -1.9222696
79 80 81 82 83 84
5.1777304 -5.2079839 11.1777304 -0.8222696 3.0920161 9.2377743
85
-7.4209044
> postscript(file="/var/www/html/rcomp/tmp/6r5nt1227728112.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 = 85
Frequency = 1
lag(myerror, k = 1) myerror
0 5.5107509 NA
1 11.6791931 5.5107509
2 8.2649074 11.6791931
3 17.4649074 8.2649074
4 1.6649074 17.4649074
5 7.5220502 1.6649074
6 2.9220502 7.5220502
7 -2.1636641 2.9220502
8 3.8220502 -2.1636641
9 8.4220502 3.8220502
10 -0.8636641 8.4220502
11 9.7820941 -0.8636641
12 14.0234154 9.7820941
13 -4.7081424 14.0234154
14 18.2775719 -4.7081424
15 -8.9224281 18.2775719
16 -11.4224281 -8.9224281
17 -1.7652852 -11.4224281
18 -0.7652852 -1.7652852
19 -10.2509995 -0.7652852
20 -2.0652852 -10.2509995
21 -0.2652852 -2.0652852
22 -11.2509995 -0.2652852
23 0.3947587 -11.2509995
24 -10.0639200 0.3947587
25 -8.2954778 -10.0639200
26 -11.1097635 -8.2954778
27 -8.1097635 -11.1097635
28 2.0902365 -8.1097635
29 -8.3526207 2.0902365
30 -8.6526207 -8.3526207
31 1.9616650 -8.6526207
32 0.7473793 1.9616650
33 -11.9526207 0.7473793
34 -4.5383350 -11.9526207
35 -9.1925768 -4.5383350
36 0.6487445 -9.1925768
37 0.8171867 0.6487445
38 -1.9970990 0.8171867
39 -5.3970990 -1.9970990
40 2.1029010 -5.3970990
41 -2.7399561 2.1029010
42 -2.5399561 -2.7399561
43 1.6743296 -2.5399561
44 -3.4399561 1.6743296
45 -12.3399561 -3.4399561
46 -1.7256704 -12.3399561
47 -5.2799122 -1.7256704
48 1.4614091 -5.2799122
49 1.5298513 1.4614091
50 -14.7844344 1.5298513
51 -1.8844344 -14.7844344
52 1.5155656 -1.8844344
53 -2.5272916 1.5155656
54 0.2727084 -2.5272916
55 5.0869941 0.2727084
56 -12.7272916 5.0869941
57 9.1727084 -12.7272916
58 1.9869941 9.1727084
59 -4.5672477 1.9869941
60 1.1740736 -4.5672477
61 3.0425158 1.1740736
62 1.8282301 3.0425158
63 8.3282301 1.8282301
64 5.6282301 8.3282301
65 9.7853730 5.6282301
66 3.5853730 9.7853730
67 8.8996587 3.5853730
68 2.4853730 8.8996587
69 7.7853730 2.4853730
70 13.2996587 7.7853730
71 -0.3748903 13.2996587
72 -5.3335690 -0.3748903
73 -4.0651268 -5.3335690
74 -0.4794125 -4.0651268
75 -1.4794125 -0.4794125
76 -1.5794125 -1.4794125
77 -1.9222696 -1.5794125
78 5.1777304 -1.9222696
79 -5.2079839 5.1777304
80 11.1777304 -5.2079839
81 -0.8222696 11.1777304
82 3.0920161 -0.8222696
83 9.2377743 3.0920161
84 -7.4209044 9.2377743
85 NA -7.4209044
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 11.6791931 5.5107509
[2,] 8.2649074 11.6791931
[3,] 17.4649074 8.2649074
[4,] 1.6649074 17.4649074
[5,] 7.5220502 1.6649074
[6,] 2.9220502 7.5220502
[7,] -2.1636641 2.9220502
[8,] 3.8220502 -2.1636641
[9,] 8.4220502 3.8220502
[10,] -0.8636641 8.4220502
[11,] 9.7820941 -0.8636641
[12,] 14.0234154 9.7820941
[13,] -4.7081424 14.0234154
[14,] 18.2775719 -4.7081424
[15,] -8.9224281 18.2775719
[16,] -11.4224281 -8.9224281
[17,] -1.7652852 -11.4224281
[18,] -0.7652852 -1.7652852
[19,] -10.2509995 -0.7652852
[20,] -2.0652852 -10.2509995
[21,] -0.2652852 -2.0652852
[22,] -11.2509995 -0.2652852
[23,] 0.3947587 -11.2509995
[24,] -10.0639200 0.3947587
[25,] -8.2954778 -10.0639200
[26,] -11.1097635 -8.2954778
[27,] -8.1097635 -11.1097635
[28,] 2.0902365 -8.1097635
[29,] -8.3526207 2.0902365
[30,] -8.6526207 -8.3526207
[31,] 1.9616650 -8.6526207
[32,] 0.7473793 1.9616650
[33,] -11.9526207 0.7473793
[34,] -4.5383350 -11.9526207
[35,] -9.1925768 -4.5383350
[36,] 0.6487445 -9.1925768
[37,] 0.8171867 0.6487445
[38,] -1.9970990 0.8171867
[39,] -5.3970990 -1.9970990
[40,] 2.1029010 -5.3970990
[41,] -2.7399561 2.1029010
[42,] -2.5399561 -2.7399561
[43,] 1.6743296 -2.5399561
[44,] -3.4399561 1.6743296
[45,] -12.3399561 -3.4399561
[46,] -1.7256704 -12.3399561
[47,] -5.2799122 -1.7256704
[48,] 1.4614091 -5.2799122
[49,] 1.5298513 1.4614091
[50,] -14.7844344 1.5298513
[51,] -1.8844344 -14.7844344
[52,] 1.5155656 -1.8844344
[53,] -2.5272916 1.5155656
[54,] 0.2727084 -2.5272916
[55,] 5.0869941 0.2727084
[56,] -12.7272916 5.0869941
[57,] 9.1727084 -12.7272916
[58,] 1.9869941 9.1727084
[59,] -4.5672477 1.9869941
[60,] 1.1740736 -4.5672477
[61,] 3.0425158 1.1740736
[62,] 1.8282301 3.0425158
[63,] 8.3282301 1.8282301
[64,] 5.6282301 8.3282301
[65,] 9.7853730 5.6282301
[66,] 3.5853730 9.7853730
[67,] 8.8996587 3.5853730
[68,] 2.4853730 8.8996587
[69,] 7.7853730 2.4853730
[70,] 13.2996587 7.7853730
[71,] -0.3748903 13.2996587
[72,] -5.3335690 -0.3748903
[73,] -4.0651268 -5.3335690
[74,] -0.4794125 -4.0651268
[75,] -1.4794125 -0.4794125
[76,] -1.5794125 -1.4794125
[77,] -1.9222696 -1.5794125
[78,] 5.1777304 -1.9222696
[79,] -5.2079839 5.1777304
[80,] 11.1777304 -5.2079839
[81,] -0.8222696 11.1777304
[82,] 3.0920161 -0.8222696
[83,] 9.2377743 3.0920161
[84,] -7.4209044 9.2377743
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 11.6791931 5.5107509
2 8.2649074 11.6791931
3 17.4649074 8.2649074
4 1.6649074 17.4649074
5 7.5220502 1.6649074
6 2.9220502 7.5220502
7 -2.1636641 2.9220502
8 3.8220502 -2.1636641
9 8.4220502 3.8220502
10 -0.8636641 8.4220502
11 9.7820941 -0.8636641
12 14.0234154 9.7820941
13 -4.7081424 14.0234154
14 18.2775719 -4.7081424
15 -8.9224281 18.2775719
16 -11.4224281 -8.9224281
17 -1.7652852 -11.4224281
18 -0.7652852 -1.7652852
19 -10.2509995 -0.7652852
20 -2.0652852 -10.2509995
21 -0.2652852 -2.0652852
22 -11.2509995 -0.2652852
23 0.3947587 -11.2509995
24 -10.0639200 0.3947587
25 -8.2954778 -10.0639200
26 -11.1097635 -8.2954778
27 -8.1097635 -11.1097635
28 2.0902365 -8.1097635
29 -8.3526207 2.0902365
30 -8.6526207 -8.3526207
31 1.9616650 -8.6526207
32 0.7473793 1.9616650
33 -11.9526207 0.7473793
34 -4.5383350 -11.9526207
35 -9.1925768 -4.5383350
36 0.6487445 -9.1925768
37 0.8171867 0.6487445
38 -1.9970990 0.8171867
39 -5.3970990 -1.9970990
40 2.1029010 -5.3970990
41 -2.7399561 2.1029010
42 -2.5399561 -2.7399561
43 1.6743296 -2.5399561
44 -3.4399561 1.6743296
45 -12.3399561 -3.4399561
46 -1.7256704 -12.3399561
47 -5.2799122 -1.7256704
48 1.4614091 -5.2799122
49 1.5298513 1.4614091
50 -14.7844344 1.5298513
51 -1.8844344 -14.7844344
52 1.5155656 -1.8844344
53 -2.5272916 1.5155656
54 0.2727084 -2.5272916
55 5.0869941 0.2727084
56 -12.7272916 5.0869941
57 9.1727084 -12.7272916
58 1.9869941 9.1727084
59 -4.5672477 1.9869941
60 1.1740736 -4.5672477
61 3.0425158 1.1740736
62 1.8282301 3.0425158
63 8.3282301 1.8282301
64 5.6282301 8.3282301
65 9.7853730 5.6282301
66 3.5853730 9.7853730
67 8.8996587 3.5853730
68 2.4853730 8.8996587
69 7.7853730 2.4853730
70 13.2996587 7.7853730
71 -0.3748903 13.2996587
72 -5.3335690 -0.3748903
73 -4.0651268 -5.3335690
74 -0.4794125 -4.0651268
75 -1.4794125 -0.4794125
76 -1.5794125 -1.4794125
77 -1.9222696 -1.5794125
78 5.1777304 -1.9222696
79 -5.2079839 5.1777304
80 11.1777304 -5.2079839
81 -0.8222696 11.1777304
82 3.0920161 -0.8222696
83 9.2377743 3.0920161
84 -7.4209044 9.2377743
> 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/73w3n1227728112.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/8gmog1227728112.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/9p4w41227728112.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/10roho1227728112.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/11gkf51227728112.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/122icp1227728112.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/13x47a1227728112.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/14r8ir1227728113.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/152i8b1227728113.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/16j91p1227728113.tab")
+ }
>
> system("convert tmp/1w2a31227728112.ps tmp/1w2a31227728112.png")
> system("convert tmp/2gm5g1227728112.ps tmp/2gm5g1227728112.png")
> system("convert tmp/3n3ho1227728112.ps tmp/3n3ho1227728112.png")
> system("convert tmp/4hymi1227728112.ps tmp/4hymi1227728112.png")
> system("convert tmp/5olhb1227728112.ps tmp/5olhb1227728112.png")
> system("convert tmp/6r5nt1227728112.ps tmp/6r5nt1227728112.png")
> system("convert tmp/73w3n1227728112.ps tmp/73w3n1227728112.png")
> system("convert tmp/8gmog1227728112.ps tmp/8gmog1227728112.png")
> system("convert tmp/9p4w41227728112.ps tmp/9p4w41227728112.png")
> system("convert tmp/10roho1227728112.ps tmp/10roho1227728112.png")
>
>
> proc.time()
user system elapsed
2.753 1.594 3.387