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(4223.4,401,4627.3,394,5175.3,372,4550.7,334,4639.3,320,5498.7,334,5031.0,400,4033.3,427,4643.5,423,4873.2,395,4608.7,373,4733.5,377,3955.6,391,4590.9,398,5127.5,393,5257.3,375,5416.9,371,5813.3,364,5261.9,400,4669.2,406,5855.8,407,5274.6,397,5516.7,389,5819.5,394,5156.0,399,5377.3,401,6386.8,396,5144.0,392,6138.5,384,5567.8,370,5822.6,380,5145.5,376,5706.6,378,6078.5,376,6074.5,373,5577.6,374,5727.5,379,6067.0,376,7069.9,371,5490.0,375,5948.3,360,6177.5,338,6890.1,352,5756.2,344,6528.8,330,6792.0,334,6657.4,333,5753.7,343,5750.9,350,5968.4,341,5871.7,320,7004.9,302,6363.4,287,6694.7,304,7101.6,370,5364.0,385,6958.6,365,6503.3,333,5316.0,313,5312.7,330,4478.0,367),dim=c(2,61),dimnames=list(c('Export','Werkloosheid'),1:61))
> y <- array(NA,dim=c(2,61),dimnames=list(c('Export','Werkloosheid'),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 = '2'
> #'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
Werkloosheid Export M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 401 4223.4 1 0 0 0 0 0 0 0 0 0 0 1
2 394 4627.3 0 1 0 0 0 0 0 0 0 0 0 2
3 372 5175.3 0 0 1 0 0 0 0 0 0 0 0 3
4 334 4550.7 0 0 0 1 0 0 0 0 0 0 0 4
5 320 4639.3 0 0 0 0 1 0 0 0 0 0 0 5
6 334 5498.7 0 0 0 0 0 1 0 0 0 0 0 6
7 400 5031.0 0 0 0 0 0 0 1 0 0 0 0 7
8 427 4033.3 0 0 0 0 0 0 0 1 0 0 0 8
9 423 4643.5 0 0 0 0 0 0 0 0 1 0 0 9
10 395 4873.2 0 0 0 0 0 0 0 0 0 1 0 10
11 373 4608.7 0 0 0 0 0 0 0 0 0 0 1 11
12 377 4733.5 0 0 0 0 0 0 0 0 0 0 0 12
13 391 3955.6 1 0 0 0 0 0 0 0 0 0 0 13
14 398 4590.9 0 1 0 0 0 0 0 0 0 0 0 14
15 393 5127.5 0 0 1 0 0 0 0 0 0 0 0 15
16 375 5257.3 0 0 0 1 0 0 0 0 0 0 0 16
17 371 5416.9 0 0 0 0 1 0 0 0 0 0 0 17
18 364 5813.3 0 0 0 0 0 1 0 0 0 0 0 18
19 400 5261.9 0 0 0 0 0 0 1 0 0 0 0 19
20 406 4669.2 0 0 0 0 0 0 0 1 0 0 0 20
21 407 5855.8 0 0 0 0 0 0 0 0 1 0 0 21
22 397 5274.6 0 0 0 0 0 0 0 0 0 1 0 22
23 389 5516.7 0 0 0 0 0 0 0 0 0 0 1 23
24 394 5819.5 0 0 0 0 0 0 0 0 0 0 0 24
25 399 5156.0 1 0 0 0 0 0 0 0 0 0 0 25
26 401 5377.3 0 1 0 0 0 0 0 0 0 0 0 26
27 396 6386.8 0 0 1 0 0 0 0 0 0 0 0 27
28 392 5144.0 0 0 0 1 0 0 0 0 0 0 0 28
29 384 6138.5 0 0 0 0 1 0 0 0 0 0 0 29
30 370 5567.8 0 0 0 0 0 1 0 0 0 0 0 30
31 380 5822.6 0 0 0 0 0 0 1 0 0 0 0 31
32 376 5145.5 0 0 0 0 0 0 0 1 0 0 0 32
33 378 5706.6 0 0 0 0 0 0 0 0 1 0 0 33
34 376 6078.5 0 0 0 0 0 0 0 0 0 1 0 34
35 373 6074.5 0 0 0 0 0 0 0 0 0 0 1 35
36 374 5577.6 0 0 0 0 0 0 0 0 0 0 0 36
37 379 5727.5 1 0 0 0 0 0 0 0 0 0 0 37
38 376 6067.0 0 1 0 0 0 0 0 0 0 0 0 38
39 371 7069.9 0 0 1 0 0 0 0 0 0 0 0 39
40 375 5490.0 0 0 0 1 0 0 0 0 0 0 0 40
41 360 5948.3 0 0 0 0 1 0 0 0 0 0 0 41
42 338 6177.5 0 0 0 0 0 1 0 0 0 0 0 42
43 352 6890.1 0 0 0 0 0 0 1 0 0 0 0 43
44 344 5756.2 0 0 0 0 0 0 0 1 0 0 0 44
45 330 6528.8 0 0 0 0 0 0 0 0 1 0 0 45
46 334 6792.0 0 0 0 0 0 0 0 0 0 1 0 46
47 333 6657.4 0 0 0 0 0 0 0 0 0 0 1 47
48 343 5753.7 0 0 0 0 0 0 0 0 0 0 0 48
49 350 5750.9 1 0 0 0 0 0 0 0 0 0 0 49
50 341 5968.4 0 1 0 0 0 0 0 0 0 0 0 50
51 320 5871.7 0 0 1 0 0 0 0 0 0 0 0 51
52 302 7004.9 0 0 0 1 0 0 0 0 0 0 0 52
53 287 6363.4 0 0 0 0 1 0 0 0 0 0 0 53
54 304 6694.7 0 0 0 0 0 1 0 0 0 0 0 54
55 370 7101.6 0 0 0 0 0 0 1 0 0 0 0 55
56 385 5364.0 0 0 0 0 0 0 0 1 0 0 0 56
57 365 6958.6 0 0 0 0 0 0 0 0 1 0 0 57
58 333 6503.3 0 0 0 0 0 0 0 0 0 1 0 58
59 313 5316.0 0 0 0 0 0 0 0 0 0 0 1 59
60 330 5312.7 0 0 0 0 0 0 0 0 0 0 0 60
61 367 4478.0 1 0 0 0 0 0 0 0 0 0 0 61
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Export M1 M2 M3 M4
4.023e+02 -2.817e-04 1.225e+01 8.045e+00 -2.354e+00 -1.624e+01
M5 M6 M7 M8 M9 M10
-2.635e+01 -2.765e+01 1.180e+01 1.975e+01 1.404e+01 1.466e+00
M11 t
-8.377e+00 -1.032e+00
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-49.475 -12.340 4.266 14.353 39.723
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.023e+02 3.181e+01 12.649 < 2e-16 ***
Export -2.817e-04 6.777e-03 -0.042 0.967015
M1 1.225e+01 1.416e+01 0.865 0.391481
M2 8.045e+00 1.464e+01 0.550 0.585150
M3 -2.354e+00 1.544e+01 -0.152 0.879500
M4 -1.624e+01 1.468e+01 -1.107 0.274041
M5 -2.635e+01 1.488e+01 -1.771 0.083033 .
M6 -2.765e+01 1.526e+01 -1.812 0.076315 .
M7 1.180e+01 1.533e+01 0.770 0.445315
M8 1.975e+01 1.466e+01 1.347 0.184578
M9 1.404e+01 1.503e+01 0.934 0.354934
M10 1.466e+00 1.492e+01 0.098 0.922110
M11 -8.377e+00 1.456e+01 -0.575 0.567858
t -1.032e+00 2.655e-01 -3.888 0.000316 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 22.9 on 47 degrees of freedom
Multiple R-squared: 0.5736, Adjusted R-squared: 0.4557
F-statistic: 4.864 on 13 and 47 DF, p-value: 2.723e-05
> 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.05514875 0.11029750 0.9448512
[2,] 0.01747354 0.03494708 0.9825265
[3,] 0.03861143 0.07722287 0.9613886
[4,] 0.41397395 0.82794790 0.5860260
[5,] 0.53779329 0.92441342 0.4622067
[6,] 0.42576894 0.85153789 0.5742311
[7,] 0.31506470 0.63012940 0.6849353
[8,] 0.22210644 0.44421288 0.7778936
[9,] 0.17020180 0.34040360 0.8297982
[10,] 0.11503872 0.23007744 0.8849613
[11,] 0.07618377 0.15236754 0.9238162
[12,] 0.07211651 0.14423302 0.9278835
[13,] 0.07856079 0.15712157 0.9214392
[14,] 0.05008301 0.10016602 0.9499170
[15,] 0.08072438 0.16144876 0.9192756
[16,] 0.19205022 0.38410044 0.8079498
[17,] 0.24421571 0.48843142 0.7557843
[18,] 0.20272179 0.40544359 0.7972782
[19,] 0.16002854 0.32005707 0.8399715
[20,] 0.11448450 0.22896901 0.8855155
[21,] 0.08205866 0.16411731 0.9179413
[22,] 0.06246069 0.12492138 0.9375393
[23,] 0.09831060 0.19662120 0.9016894
[24,] 0.09215900 0.18431800 0.9078410
[25,] 0.33243840 0.66487679 0.6675616
[26,] 0.62582266 0.74835469 0.3741773
[27,] 0.52283977 0.95432046 0.4771602
[28,] 0.64253519 0.71492962 0.3574648
> postscript(file="/var/www/html/rcomp/tmp/1memd1262191055.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/2hjl41262191055.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/377kl1262191055.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/49knp1262191055.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/5kan31262191055.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
-11.3222819 -12.9729741 -23.3876415 -46.6405354 -49.4752740 -32.9033335
7 8 9 10 11 12
-5.4551183 14.3533621 17.2590298 2.9334205 -8.2651257 -11.5749519
13 14 15 16 17 18
-9.0096963 3.4048062 9.9869270 6.9465780 14.1318430 9.5733377
19 20 21 22 23 24
6.9979712 5.9205565 13.9886193 17.4345468 20.3787302 18.1190538
25 26 27 28 29 30
11.7165405 19.0144026 25.7297583 36.3026925 39.7231825 27.8922061
31 32 33 34 35 36
-0.4560213 -11.5572149 -2.6653807 9.0490735 16.9239207 10.4389365
37 38 39 40 41 42
4.2655908 6.5967546 13.3102509 31.7882104 28.0576312 8.4520189
43 44 45 46 47 48
-15.7672277 -30.9971204 -38.0456980 -20.3618690 -10.5238172 -8.1234133
49 50 51 52 53 54
-12.3397807 -16.0429893 -25.6392947 -28.3969454 -32.4373827 -13.0142293
55 56 57 58 59 60
14.6803961 22.2804168 9.4634297 -9.0551718 -18.5137080 -8.8596251
61
16.6896276
> postscript(file="/var/www/html/rcomp/tmp/6l7lo1262191055.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 -11.3222819 NA
1 -12.9729741 -11.3222819
2 -23.3876415 -12.9729741
3 -46.6405354 -23.3876415
4 -49.4752740 -46.6405354
5 -32.9033335 -49.4752740
6 -5.4551183 -32.9033335
7 14.3533621 -5.4551183
8 17.2590298 14.3533621
9 2.9334205 17.2590298
10 -8.2651257 2.9334205
11 -11.5749519 -8.2651257
12 -9.0096963 -11.5749519
13 3.4048062 -9.0096963
14 9.9869270 3.4048062
15 6.9465780 9.9869270
16 14.1318430 6.9465780
17 9.5733377 14.1318430
18 6.9979712 9.5733377
19 5.9205565 6.9979712
20 13.9886193 5.9205565
21 17.4345468 13.9886193
22 20.3787302 17.4345468
23 18.1190538 20.3787302
24 11.7165405 18.1190538
25 19.0144026 11.7165405
26 25.7297583 19.0144026
27 36.3026925 25.7297583
28 39.7231825 36.3026925
29 27.8922061 39.7231825
30 -0.4560213 27.8922061
31 -11.5572149 -0.4560213
32 -2.6653807 -11.5572149
33 9.0490735 -2.6653807
34 16.9239207 9.0490735
35 10.4389365 16.9239207
36 4.2655908 10.4389365
37 6.5967546 4.2655908
38 13.3102509 6.5967546
39 31.7882104 13.3102509
40 28.0576312 31.7882104
41 8.4520189 28.0576312
42 -15.7672277 8.4520189
43 -30.9971204 -15.7672277
44 -38.0456980 -30.9971204
45 -20.3618690 -38.0456980
46 -10.5238172 -20.3618690
47 -8.1234133 -10.5238172
48 -12.3397807 -8.1234133
49 -16.0429893 -12.3397807
50 -25.6392947 -16.0429893
51 -28.3969454 -25.6392947
52 -32.4373827 -28.3969454
53 -13.0142293 -32.4373827
54 14.6803961 -13.0142293
55 22.2804168 14.6803961
56 9.4634297 22.2804168
57 -9.0551718 9.4634297
58 -18.5137080 -9.0551718
59 -8.8596251 -18.5137080
60 16.6896276 -8.8596251
61 NA 16.6896276
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -12.9729741 -11.3222819
[2,] -23.3876415 -12.9729741
[3,] -46.6405354 -23.3876415
[4,] -49.4752740 -46.6405354
[5,] -32.9033335 -49.4752740
[6,] -5.4551183 -32.9033335
[7,] 14.3533621 -5.4551183
[8,] 17.2590298 14.3533621
[9,] 2.9334205 17.2590298
[10,] -8.2651257 2.9334205
[11,] -11.5749519 -8.2651257
[12,] -9.0096963 -11.5749519
[13,] 3.4048062 -9.0096963
[14,] 9.9869270 3.4048062
[15,] 6.9465780 9.9869270
[16,] 14.1318430 6.9465780
[17,] 9.5733377 14.1318430
[18,] 6.9979712 9.5733377
[19,] 5.9205565 6.9979712
[20,] 13.9886193 5.9205565
[21,] 17.4345468 13.9886193
[22,] 20.3787302 17.4345468
[23,] 18.1190538 20.3787302
[24,] 11.7165405 18.1190538
[25,] 19.0144026 11.7165405
[26,] 25.7297583 19.0144026
[27,] 36.3026925 25.7297583
[28,] 39.7231825 36.3026925
[29,] 27.8922061 39.7231825
[30,] -0.4560213 27.8922061
[31,] -11.5572149 -0.4560213
[32,] -2.6653807 -11.5572149
[33,] 9.0490735 -2.6653807
[34,] 16.9239207 9.0490735
[35,] 10.4389365 16.9239207
[36,] 4.2655908 10.4389365
[37,] 6.5967546 4.2655908
[38,] 13.3102509 6.5967546
[39,] 31.7882104 13.3102509
[40,] 28.0576312 31.7882104
[41,] 8.4520189 28.0576312
[42,] -15.7672277 8.4520189
[43,] -30.9971204 -15.7672277
[44,] -38.0456980 -30.9971204
[45,] -20.3618690 -38.0456980
[46,] -10.5238172 -20.3618690
[47,] -8.1234133 -10.5238172
[48,] -12.3397807 -8.1234133
[49,] -16.0429893 -12.3397807
[50,] -25.6392947 -16.0429893
[51,] -28.3969454 -25.6392947
[52,] -32.4373827 -28.3969454
[53,] -13.0142293 -32.4373827
[54,] 14.6803961 -13.0142293
[55,] 22.2804168 14.6803961
[56,] 9.4634297 22.2804168
[57,] -9.0551718 9.4634297
[58,] -18.5137080 -9.0551718
[59,] -8.8596251 -18.5137080
[60,] 16.6896276 -8.8596251
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -12.9729741 -11.3222819
2 -23.3876415 -12.9729741
3 -46.6405354 -23.3876415
4 -49.4752740 -46.6405354
5 -32.9033335 -49.4752740
6 -5.4551183 -32.9033335
7 14.3533621 -5.4551183
8 17.2590298 14.3533621
9 2.9334205 17.2590298
10 -8.2651257 2.9334205
11 -11.5749519 -8.2651257
12 -9.0096963 -11.5749519
13 3.4048062 -9.0096963
14 9.9869270 3.4048062
15 6.9465780 9.9869270
16 14.1318430 6.9465780
17 9.5733377 14.1318430
18 6.9979712 9.5733377
19 5.9205565 6.9979712
20 13.9886193 5.9205565
21 17.4345468 13.9886193
22 20.3787302 17.4345468
23 18.1190538 20.3787302
24 11.7165405 18.1190538
25 19.0144026 11.7165405
26 25.7297583 19.0144026
27 36.3026925 25.7297583
28 39.7231825 36.3026925
29 27.8922061 39.7231825
30 -0.4560213 27.8922061
31 -11.5572149 -0.4560213
32 -2.6653807 -11.5572149
33 9.0490735 -2.6653807
34 16.9239207 9.0490735
35 10.4389365 16.9239207
36 4.2655908 10.4389365
37 6.5967546 4.2655908
38 13.3102509 6.5967546
39 31.7882104 13.3102509
40 28.0576312 31.7882104
41 8.4520189 28.0576312
42 -15.7672277 8.4520189
43 -30.9971204 -15.7672277
44 -38.0456980 -30.9971204
45 -20.3618690 -38.0456980
46 -10.5238172 -20.3618690
47 -8.1234133 -10.5238172
48 -12.3397807 -8.1234133
49 -16.0429893 -12.3397807
50 -25.6392947 -16.0429893
51 -28.3969454 -25.6392947
52 -32.4373827 -28.3969454
53 -13.0142293 -32.4373827
54 14.6803961 -13.0142293
55 22.2804168 14.6803961
56 9.4634297 22.2804168
57 -9.0551718 9.4634297
58 -18.5137080 -9.0551718
59 -8.8596251 -18.5137080
60 16.6896276 -8.8596251
> 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/746by1262191055.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/8ppra1262191055.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/9gl5t1262191055.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/10x9tz1262191055.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/1100rt1262191055.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/12m2od1262191055.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/13uxw81262191056.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/14ljau1262191056.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/15zkxu1262191056.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/16hq071262191056.tab")
+ }
>
> try(system("convert tmp/1memd1262191055.ps tmp/1memd1262191055.png",intern=TRUE))
character(0)
> try(system("convert tmp/2hjl41262191055.ps tmp/2hjl41262191055.png",intern=TRUE))
character(0)
> try(system("convert tmp/377kl1262191055.ps tmp/377kl1262191055.png",intern=TRUE))
character(0)
> try(system("convert tmp/49knp1262191055.ps tmp/49knp1262191055.png",intern=TRUE))
character(0)
> try(system("convert tmp/5kan31262191055.ps tmp/5kan31262191055.png",intern=TRUE))
character(0)
> try(system("convert tmp/6l7lo1262191055.ps tmp/6l7lo1262191055.png",intern=TRUE))
character(0)
> try(system("convert tmp/746by1262191055.ps tmp/746by1262191055.png",intern=TRUE))
character(0)
> try(system("convert tmp/8ppra1262191055.ps tmp/8ppra1262191055.png",intern=TRUE))
character(0)
> try(system("convert tmp/9gl5t1262191055.ps tmp/9gl5t1262191055.png",intern=TRUE))
character(0)
> try(system("convert tmp/10x9tz1262191055.ps tmp/10x9tz1262191055.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
2.434 1.576 3.207