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(118.7,0,110.1,0,110.3,0,112.9,0,102.2,0,99.4,0,116.1,0,103.8,0,101.8,0,113.7,0,89.7,0,99.5,0,122.9,0,108.6,0,114.4,0,110.5,0,104.1,0,103.6,0,121.6,0,101.1,0,116.0,0,120.1,0,96.0,0,105.0,0,124.7,0,123.9,0,123.6,0,114.8,0,108.8,0,106.1,0,123.2,0,106.2,0,115.2,0,120.6,0,109.5,0,114.4,0,121.4,0,129.5,0,124.3,0,112.6,0,125.1,1,117.9,1,116.4,1,126.4,1,93.3,1,102.9,1,97.8,1,97.1,1,110.7,1,109.3,1,103.2,1,106.2,1,81.3,1,84.5,1,92.7,1,85.0,1,79.1,1,92.6,1,78.1,1,76.9,1,92.5,1),dim=c(2,61),dimnames=list(c('Y(t)_Bruto_index_consumptiegoederen','Dummyvariabele'),1:61))
> y <- array(NA,dim=c(2,61),dimnames=list(c('Y(t)_Bruto_index_consumptiegoederen','Dummyvariabele'),1:61))
> 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(t)_Bruto_index_consumptiegoederen Dummyvariabele M1 M2 M3 M4 M5 M6 M7 M8
1 118.7 0 1 0 0 0 0 0 0 0
2 110.1 0 0 1 0 0 0 0 0 0
3 110.3 0 0 0 1 0 0 0 0 0
4 112.9 0 0 0 0 1 0 0 0 0
5 102.2 0 0 0 0 0 1 0 0 0
6 99.4 0 0 0 0 0 0 1 0 0
7 116.1 0 0 0 0 0 0 0 1 0
8 103.8 0 0 0 0 0 0 0 0 1
9 101.8 0 0 0 0 0 0 0 0 0
10 113.7 0 0 0 0 0 0 0 0 0
11 89.7 0 0 0 0 0 0 0 0 0
12 99.5 0 0 0 0 0 0 0 0 0
13 122.9 0 1 0 0 0 0 0 0 0
14 108.6 0 0 1 0 0 0 0 0 0
15 114.4 0 0 0 1 0 0 0 0 0
16 110.5 0 0 0 0 1 0 0 0 0
17 104.1 0 0 0 0 0 1 0 0 0
18 103.6 0 0 0 0 0 0 1 0 0
19 121.6 0 0 0 0 0 0 0 1 0
20 101.1 0 0 0 0 0 0 0 0 1
21 116.0 0 0 0 0 0 0 0 0 0
22 120.1 0 0 0 0 0 0 0 0 0
23 96.0 0 0 0 0 0 0 0 0 0
24 105.0 0 0 0 0 0 0 0 0 0
25 124.7 0 1 0 0 0 0 0 0 0
26 123.9 0 0 1 0 0 0 0 0 0
27 123.6 0 0 0 1 0 0 0 0 0
28 114.8 0 0 0 0 1 0 0 0 0
29 108.8 0 0 0 0 0 1 0 0 0
30 106.1 0 0 0 0 0 0 1 0 0
31 123.2 0 0 0 0 0 0 0 1 0
32 106.2 0 0 0 0 0 0 0 0 1
33 115.2 0 0 0 0 0 0 0 0 0
34 120.6 0 0 0 0 0 0 0 0 0
35 109.5 0 0 0 0 0 0 0 0 0
36 114.4 0 0 0 0 0 0 0 0 0
37 121.4 0 1 0 0 0 0 0 0 0
38 129.5 0 0 1 0 0 0 0 0 0
39 124.3 0 0 0 1 0 0 0 0 0
40 112.6 0 0 0 0 1 0 0 0 0
41 125.1 1 0 0 0 0 1 0 0 0
42 117.9 1 0 0 0 0 0 1 0 0
43 116.4 1 0 0 0 0 0 0 1 0
44 126.4 1 0 0 0 0 0 0 0 1
45 93.3 1 0 0 0 0 0 0 0 0
46 102.9 1 0 0 0 0 0 0 0 0
47 97.8 1 0 0 0 0 0 0 0 0
48 97.1 1 0 0 0 0 0 0 0 0
49 110.7 1 1 0 0 0 0 0 0 0
50 109.3 1 0 1 0 0 0 0 0 0
51 103.2 1 0 0 1 0 0 0 0 0
52 106.2 1 0 0 0 1 0 0 0 0
53 81.3 1 0 0 0 0 1 0 0 0
54 84.5 1 0 0 0 0 0 1 0 0
55 92.7 1 0 0 0 0 0 0 1 0
56 85.0 1 0 0 0 0 0 0 0 1
57 79.1 1 0 0 0 0 0 0 0 0
58 92.6 1 0 0 0 0 0 0 0 0
59 78.1 1 0 0 0 0 0 0 0 0
60 76.9 1 0 0 0 0 0 0 0 0
61 92.5 1 1 0 0 0 0 0 0 0
M9 M10 M11 t
1 0 0 0 1
2 0 0 0 2
3 0 0 0 3
4 0 0 0 4
5 0 0 0 5
6 0 0 0 6
7 0 0 0 7
8 0 0 0 8
9 1 0 0 9
10 0 1 0 10
11 0 0 1 11
12 0 0 0 12
13 0 0 0 13
14 0 0 0 14
15 0 0 0 15
16 0 0 0 16
17 0 0 0 17
18 0 0 0 18
19 0 0 0 19
20 0 0 0 20
21 1 0 0 21
22 0 1 0 22
23 0 0 1 23
24 0 0 0 24
25 0 0 0 25
26 0 0 0 26
27 0 0 0 27
28 0 0 0 28
29 0 0 0 29
30 0 0 0 30
31 0 0 0 31
32 0 0 0 32
33 1 0 0 33
34 0 1 0 34
35 0 0 1 35
36 0 0 0 36
37 0 0 0 37
38 0 0 0 38
39 0 0 0 39
40 0 0 0 40
41 0 0 0 41
42 0 0 0 42
43 0 0 0 43
44 0 0 0 44
45 1 0 0 45
46 0 1 0 46
47 0 0 1 47
48 0 0 0 48
49 0 0 0 49
50 0 0 0 50
51 0 0 0 51
52 0 0 0 52
53 0 0 0 53
54 0 0 0 54
55 0 0 0 55
56 0 0 0 56
57 1 0 0 57
58 0 1 0 58
59 0 0 1 59
60 0 0 0 60
61 0 0 0 61
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Dummyvariabele M1 M2 M3
102.03915 -13.96094 15.93444 15.49815 14.31912
M4 M5 M6 M7 M8
10.50008 6.13324 4.07420 15.71517 6.15614
M9 M10 M11 t
2.67710 11.51807 -4.30097 0.05903
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-16.0403 -6.2354 -0.4476 5.0356 29.5682
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 102.03915 5.94029 17.177 < 2e-16 ***
Dummyvariabele -13.96094 5.18115 -2.695 0.00974 **
M1 15.93444 6.49657 2.453 0.01794 *
M2 15.49815 6.81552 2.274 0.02758 *
M3 14.31912 6.80585 2.104 0.04076 *
M4 10.50008 6.79903 1.544 0.12921
M5 6.13324 6.83882 0.897 0.37438
M6 4.07420 6.82027 0.597 0.55313
M7 15.71517 6.80453 2.310 0.02536 *
M8 6.15614 6.79162 0.906 0.36933
M9 2.67710 6.78157 0.395 0.69480
M10 11.51807 6.77438 1.700 0.09570 .
M11 -4.30097 6.77006 -0.635 0.52832
t 0.05903 0.13962 0.423 0.67436
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 10.7 on 47 degrees of freedom
Multiple R-squared: 0.4843, Adjusted R-squared: 0.3416
F-statistic: 3.395 on 13 and 47 DF, p-value: 0.001019
> 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.0169415752 0.0338831504 0.9830584
[2,] 0.0049352355 0.0098704710 0.9950648
[3,] 0.0016637691 0.0033275382 0.9983362
[4,] 0.0012390461 0.0024780922 0.9987610
[5,] 0.0069966989 0.0139933979 0.9930033
[6,] 0.0028622566 0.0057245132 0.9971377
[7,] 0.0019019998 0.0038039995 0.9980980
[8,] 0.0010160400 0.0020320801 0.9989840
[9,] 0.0004462129 0.0008924259 0.9995538
[10,] 0.0012886369 0.0025772738 0.9987114
[11,] 0.0009864987 0.0019729973 0.9990135
[12,] 0.0014060697 0.0028121393 0.9985939
[13,] 0.0009953733 0.0019907465 0.9990046
[14,] 0.0010338198 0.0020676397 0.9989662
[15,] 0.0004722714 0.0009445429 0.9995277
[16,] 0.0027841989 0.0055683978 0.9972158
[17,] 0.0013411072 0.0026822144 0.9986589
[18,] 0.0006495064 0.0012990129 0.9993505
[19,] 0.0019920884 0.0039841769 0.9980079
[20,] 0.0013231375 0.0026462750 0.9986769
[21,] 0.0029101959 0.0058203919 0.9970898
[22,] 0.0025235663 0.0050471327 0.9974764
[23,] 0.0019282733 0.0038565465 0.9980717
[24,] 0.0018333300 0.0036666601 0.9981667
[25,] 0.0092949116 0.0185898233 0.9907051
[26,] 0.0123217299 0.0246434598 0.9876783
[27,] 0.0315765930 0.0631531859 0.9684234
[28,] 0.7164891010 0.5670217980 0.2835109
> postscript(file="/var/www/html/rcomp/tmp/1rlgw1260543192.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/25ht61260543192.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/3aa531260543192.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/40q3f1260543192.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/5yqp21260543192.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 = 61
Frequency = 1
1 2 3 4 5 6
0.6673713 -7.5553738 -6.2353738 0.1246262 -6.2675613 -7.0675613
7 8 9 10 11 12
-2.0675613 -4.8675613 -3.4475613 -0.4475613 -8.6875613 -3.2475613
13 14 15 16 17 18
4.1589645 -9.7637806 -2.8437806 -2.9837806 -5.0759681 -3.5759681
19 20 21 22 23 24
2.7240319 -8.2759681 10.0440319 5.2440319 -3.0959681 1.5440319
25 26 27 28 29 30
5.2505576 4.8278125 5.6478125 0.6078125 -1.0843750 -1.7843750
31 32 33 34 35 36
3.6156250 -3.8843750 8.5356250 5.0356250 9.6956250 10.2356250
37 38 39 40 41 42
1.2421507 9.7194056 5.6394056 -2.3005944 28.4681556 23.2681556
43 44 45 46 47 48
10.0681556 29.5681556 -0.1118444 0.5881556 11.2481556 6.1881556
49 50 51 52 53 54
3.7946814 2.7719363 -2.2080637 4.5519363 -16.0402512 -10.8402512
55 56 57 58 59 60
-14.3402512 -12.5402512 -15.0202512 -10.4202512 -9.1602512 -14.7202512
61
-15.1137255
> postscript(file="/var/www/html/rcomp/tmp/6d5w91260543192.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 = 61
Frequency = 1
lag(myerror, k = 1) myerror
0 0.6673713 NA
1 -7.5553738 0.6673713
2 -6.2353738 -7.5553738
3 0.1246262 -6.2353738
4 -6.2675613 0.1246262
5 -7.0675613 -6.2675613
6 -2.0675613 -7.0675613
7 -4.8675613 -2.0675613
8 -3.4475613 -4.8675613
9 -0.4475613 -3.4475613
10 -8.6875613 -0.4475613
11 -3.2475613 -8.6875613
12 4.1589645 -3.2475613
13 -9.7637806 4.1589645
14 -2.8437806 -9.7637806
15 -2.9837806 -2.8437806
16 -5.0759681 -2.9837806
17 -3.5759681 -5.0759681
18 2.7240319 -3.5759681
19 -8.2759681 2.7240319
20 10.0440319 -8.2759681
21 5.2440319 10.0440319
22 -3.0959681 5.2440319
23 1.5440319 -3.0959681
24 5.2505576 1.5440319
25 4.8278125 5.2505576
26 5.6478125 4.8278125
27 0.6078125 5.6478125
28 -1.0843750 0.6078125
29 -1.7843750 -1.0843750
30 3.6156250 -1.7843750
31 -3.8843750 3.6156250
32 8.5356250 -3.8843750
33 5.0356250 8.5356250
34 9.6956250 5.0356250
35 10.2356250 9.6956250
36 1.2421507 10.2356250
37 9.7194056 1.2421507
38 5.6394056 9.7194056
39 -2.3005944 5.6394056
40 28.4681556 -2.3005944
41 23.2681556 28.4681556
42 10.0681556 23.2681556
43 29.5681556 10.0681556
44 -0.1118444 29.5681556
45 0.5881556 -0.1118444
46 11.2481556 0.5881556
47 6.1881556 11.2481556
48 3.7946814 6.1881556
49 2.7719363 3.7946814
50 -2.2080637 2.7719363
51 4.5519363 -2.2080637
52 -16.0402512 4.5519363
53 -10.8402512 -16.0402512
54 -14.3402512 -10.8402512
55 -12.5402512 -14.3402512
56 -15.0202512 -12.5402512
57 -10.4202512 -15.0202512
58 -9.1602512 -10.4202512
59 -14.7202512 -9.1602512
60 -15.1137255 -14.7202512
61 NA -15.1137255
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -7.5553738 0.6673713
[2,] -6.2353738 -7.5553738
[3,] 0.1246262 -6.2353738
[4,] -6.2675613 0.1246262
[5,] -7.0675613 -6.2675613
[6,] -2.0675613 -7.0675613
[7,] -4.8675613 -2.0675613
[8,] -3.4475613 -4.8675613
[9,] -0.4475613 -3.4475613
[10,] -8.6875613 -0.4475613
[11,] -3.2475613 -8.6875613
[12,] 4.1589645 -3.2475613
[13,] -9.7637806 4.1589645
[14,] -2.8437806 -9.7637806
[15,] -2.9837806 -2.8437806
[16,] -5.0759681 -2.9837806
[17,] -3.5759681 -5.0759681
[18,] 2.7240319 -3.5759681
[19,] -8.2759681 2.7240319
[20,] 10.0440319 -8.2759681
[21,] 5.2440319 10.0440319
[22,] -3.0959681 5.2440319
[23,] 1.5440319 -3.0959681
[24,] 5.2505576 1.5440319
[25,] 4.8278125 5.2505576
[26,] 5.6478125 4.8278125
[27,] 0.6078125 5.6478125
[28,] -1.0843750 0.6078125
[29,] -1.7843750 -1.0843750
[30,] 3.6156250 -1.7843750
[31,] -3.8843750 3.6156250
[32,] 8.5356250 -3.8843750
[33,] 5.0356250 8.5356250
[34,] 9.6956250 5.0356250
[35,] 10.2356250 9.6956250
[36,] 1.2421507 10.2356250
[37,] 9.7194056 1.2421507
[38,] 5.6394056 9.7194056
[39,] -2.3005944 5.6394056
[40,] 28.4681556 -2.3005944
[41,] 23.2681556 28.4681556
[42,] 10.0681556 23.2681556
[43,] 29.5681556 10.0681556
[44,] -0.1118444 29.5681556
[45,] 0.5881556 -0.1118444
[46,] 11.2481556 0.5881556
[47,] 6.1881556 11.2481556
[48,] 3.7946814 6.1881556
[49,] 2.7719363 3.7946814
[50,] -2.2080637 2.7719363
[51,] 4.5519363 -2.2080637
[52,] -16.0402512 4.5519363
[53,] -10.8402512 -16.0402512
[54,] -14.3402512 -10.8402512
[55,] -12.5402512 -14.3402512
[56,] -15.0202512 -12.5402512
[57,] -10.4202512 -15.0202512
[58,] -9.1602512 -10.4202512
[59,] -14.7202512 -9.1602512
[60,] -15.1137255 -14.7202512
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -7.5553738 0.6673713
2 -6.2353738 -7.5553738
3 0.1246262 -6.2353738
4 -6.2675613 0.1246262
5 -7.0675613 -6.2675613
6 -2.0675613 -7.0675613
7 -4.8675613 -2.0675613
8 -3.4475613 -4.8675613
9 -0.4475613 -3.4475613
10 -8.6875613 -0.4475613
11 -3.2475613 -8.6875613
12 4.1589645 -3.2475613
13 -9.7637806 4.1589645
14 -2.8437806 -9.7637806
15 -2.9837806 -2.8437806
16 -5.0759681 -2.9837806
17 -3.5759681 -5.0759681
18 2.7240319 -3.5759681
19 -8.2759681 2.7240319
20 10.0440319 -8.2759681
21 5.2440319 10.0440319
22 -3.0959681 5.2440319
23 1.5440319 -3.0959681
24 5.2505576 1.5440319
25 4.8278125 5.2505576
26 5.6478125 4.8278125
27 0.6078125 5.6478125
28 -1.0843750 0.6078125
29 -1.7843750 -1.0843750
30 3.6156250 -1.7843750
31 -3.8843750 3.6156250
32 8.5356250 -3.8843750
33 5.0356250 8.5356250
34 9.6956250 5.0356250
35 10.2356250 9.6956250
36 1.2421507 10.2356250
37 9.7194056 1.2421507
38 5.6394056 9.7194056
39 -2.3005944 5.6394056
40 28.4681556 -2.3005944
41 23.2681556 28.4681556
42 10.0681556 23.2681556
43 29.5681556 10.0681556
44 -0.1118444 29.5681556
45 0.5881556 -0.1118444
46 11.2481556 0.5881556
47 6.1881556 11.2481556
48 3.7946814 6.1881556
49 2.7719363 3.7946814
50 -2.2080637 2.7719363
51 4.5519363 -2.2080637
52 -16.0402512 4.5519363
53 -10.8402512 -16.0402512
54 -14.3402512 -10.8402512
55 -12.5402512 -14.3402512
56 -15.0202512 -12.5402512
57 -10.4202512 -15.0202512
58 -9.1602512 -10.4202512
59 -14.7202512 -9.1602512
60 -15.1137255 -14.7202512
> 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/78e1j1260543192.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/8kmfl1260543192.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/9hgko1260543192.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/1025sa1260543192.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/11tmq31260543192.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/12hhb71260543192.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/13bfwm1260543192.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/149doh1260543192.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/15mmhn1260543192.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/16v0jh1260543192.tab")
+ }
>
> system("convert tmp/1rlgw1260543192.ps tmp/1rlgw1260543192.png")
> system("convert tmp/25ht61260543192.ps tmp/25ht61260543192.png")
> system("convert tmp/3aa531260543192.ps tmp/3aa531260543192.png")
> system("convert tmp/40q3f1260543192.ps tmp/40q3f1260543192.png")
> system("convert tmp/5yqp21260543192.ps tmp/5yqp21260543192.png")
> system("convert tmp/6d5w91260543192.ps tmp/6d5w91260543192.png")
> system("convert tmp/78e1j1260543192.ps tmp/78e1j1260543192.png")
> system("convert tmp/8kmfl1260543192.ps tmp/8kmfl1260543192.png")
> system("convert tmp/9hgko1260543192.ps tmp/9hgko1260543192.png")
> system("convert tmp/1025sa1260543192.ps tmp/1025sa1260543192.png")
>
>
> proc.time()
user system elapsed
2.470 1.593 8.233