R version 2.15.2 (2012-10-26) -- "Trick or Treat"
Copyright (C) 2012 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: i686-pc-linux-gnu (32-bit)
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(426
+ ,7.1
+ ,3.2
+ ,24776
+ ,3
+ ,396
+ ,7.2
+ ,2.9
+ ,19814
+ ,3
+ ,458
+ ,7.2
+ ,2.7
+ ,12738
+ ,3
+ ,315
+ ,7.1
+ ,3.1
+ ,31566
+ ,3
+ ,337
+ ,6.9
+ ,2.7
+ ,30111
+ ,3
+ ,386
+ ,6.8
+ ,2.6
+ ,30019
+ ,3
+ ,352
+ ,6.8
+ ,1.8
+ ,31934
+ ,3
+ ,384
+ ,6.8
+ ,2.3
+ ,25826
+ ,3
+ ,439
+ ,6.9
+ ,2.2
+ ,26835
+ ,3.18
+ ,397
+ ,7.1
+ ,1.8
+ ,20205
+ ,3.25
+ ,453
+ ,7.2
+ ,1.4
+ ,17789
+ ,3.25
+ ,364
+ ,7.2
+ ,0.3
+ ,20520
+ ,3.23
+ ,367
+ ,7.1
+ ,0.8
+ ,22518
+ ,2.92
+ ,474
+ ,7.1
+ ,-0.5
+ ,15572
+ ,2.25
+ ,373
+ ,7.2
+ ,-2.2
+ ,11509
+ ,1.62
+ ,404
+ ,7.5
+ ,-2.9
+ ,25447
+ ,1
+ ,385
+ ,7.7
+ ,-5.1
+ ,24090
+ ,0.66
+ ,365
+ ,7.8
+ ,-7.2
+ ,27786
+ ,0.31
+ ,366
+ ,7.7
+ ,-7.9
+ ,26195
+ ,0.25
+ ,421
+ ,7.7
+ ,-10.9
+ ,20516
+ ,0.25
+ ,354
+ ,7.8
+ ,-12.7
+ ,22759
+ ,0.25
+ ,367
+ ,8
+ ,-14
+ ,19028
+ ,0.25
+ ,413
+ ,8.1
+ ,-15.6
+ ,16971
+ ,0.25
+ ,362
+ ,8.1
+ ,-16
+ ,20036
+ ,0.25
+ ,385
+ ,8
+ ,-17.2
+ ,22485
+ ,0.25
+ ,343
+ ,8.1
+ ,-17.6
+ ,18730
+ ,0.25
+ ,369
+ ,8.2
+ ,-15.5
+ ,14538
+ ,0.25
+ ,363
+ ,8.4
+ ,-13.7
+ ,27561
+ ,0.25
+ ,318
+ ,8.5
+ ,-11.4
+ ,25985
+ ,0.25
+ ,393
+ ,8.5
+ ,-9.2
+ ,34670
+ ,0.25
+ ,325
+ ,8.5
+ ,-6.3
+ ,32066
+ ,0.25
+ ,403
+ ,8.5
+ ,-3.1
+ ,27186
+ ,0.25
+ ,392
+ ,8.5
+ ,0
+ ,29586
+ ,0.25
+ ,409
+ ,8.4
+ ,3
+ ,21359
+ ,0.25
+ ,485
+ ,8.3
+ ,5.4
+ ,21553
+ ,0.25
+ ,423
+ ,8.2
+ ,7.6
+ ,19573
+ ,0.25
+ ,428
+ ,8.1
+ ,9.7
+ ,24256
+ ,0.25
+ ,431
+ ,7.9
+ ,12
+ ,22380
+ ,0.25
+ ,416
+ ,7.6
+ ,11.6
+ ,16167
+ ,0.25
+ ,330
+ ,7.3
+ ,10
+ ,27297
+ ,0.25
+ ,314
+ ,7.1
+ ,10.8
+ ,28287
+ ,0.25
+ ,345
+ ,7
+ ,11.3
+ ,33474
+ ,0.39
+ ,365
+ ,7.1
+ ,10.1
+ ,28229
+ ,0.5
+ ,417
+ ,7.1
+ ,9.4
+ ,28785
+ ,0.5
+ ,356
+ ,7.1
+ ,9.6
+ ,25597
+ ,0.65
+ ,477
+ ,7.3
+ ,7.9
+ ,18130
+ ,0.75
+ ,423
+ ,7.3
+ ,7.3
+ ,20198
+ ,0.75
+ ,386
+ ,7.3
+ ,6.2
+ ,22849
+ ,0.75
+ ,390
+ ,7.2
+ ,4.9
+ ,23118
+ ,0.57
+ ,407
+ ,7.2
+ ,3.6
+ ,21925
+ ,0.36
+ ,398
+ ,7.1
+ ,2.9
+ ,20801
+ ,0.25
+ ,327
+ ,7.1
+ ,3.1
+ ,18785
+ ,0.25
+ ,304
+ ,7.1
+ ,1.7
+ ,20659
+ ,0.25
+ ,378
+ ,7.2
+ ,0.6
+ ,29367
+ ,0.25
+ ,311
+ ,7.3
+ ,-0.4
+ ,23992
+ ,0.25
+ ,376
+ ,7.4
+ ,-1.1
+ ,20645
+ ,0.25
+ ,340
+ ,7.4
+ ,-2.9
+ ,22356
+ ,0.08
+ ,383
+ ,7.5
+ ,-2.8
+ ,17902
+ ,0
+ ,467
+ ,7.4
+ ,-3
+ ,15879
+ ,0
+ ,439
+ ,7.4
+ ,-3.2
+ ,16963
+ ,0)
+ ,dim=c(5
+ ,60)
+ ,dimnames=list(c('bouwvergunningen'
+ ,'werkloosheidsgraad'
+ ,'uitvoer'
+ ,'personenwagens'
+ ,'rentetarieven')
+ ,1:60))
> y <- array(NA,dim=c(5,60),dimnames=list(c('bouwvergunningen','werkloosheidsgraad','uitvoer','personenwagens','rentetarieven'),1:60))
> for (i in 1:dim(x)[1])
+ {
+ for (j in 1:dim(x)[2])
+ {
+ y[i,j] <- as.numeric(x[i,j])
+ }
+ }
> par3 = 'No Linear Trend'
> par2 = 'Include Monthly Dummies'
> par1 = '1'
> par3 <- 'No Linear Trend'
> par2 <- 'Include Monthly Dummies'
> par1 <- '1'
> #'GNU S' R Code compiled by R2WASP v. 1.0.44 ()
> #Author: Prof. Dr. P. Wessa
> #To cite this work: AUTHOR(S), (YEAR), YOUR SOFTWARE TITLE (vNUMBER) in Free Statistics Software (v$_version), Office for Research Development and Education, URL http://www.wessa.net/rwasp_YOURPAGE.wasp/
> #Source of accompanying publication: Office for Research, Development, and Education
> #Technical description: Write here your technical program description (don't use hard returns!)
> library(lattice)
> library(lmtest)
Loading required package: zoo
Attaching package: 'zoo'
The following object(s) are masked from 'package:base':
as.Date, 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
bouwvergunningen werkloosheidsgraad uitvoer personenwagens rentetarieven M1
1 426 7.1 3.2 24776 3.00 1
2 396 7.2 2.9 19814 3.00 0
3 458 7.2 2.7 12738 3.00 0
4 315 7.1 3.1 31566 3.00 0
5 337 6.9 2.7 30111 3.00 0
6 386 6.8 2.6 30019 3.00 0
7 352 6.8 1.8 31934 3.00 0
8 384 6.8 2.3 25826 3.00 0
9 439 6.9 2.2 26835 3.18 0
10 397 7.1 1.8 20205 3.25 0
11 453 7.2 1.4 17789 3.25 0
12 364 7.2 0.3 20520 3.23 0
13 367 7.1 0.8 22518 2.92 1
14 474 7.1 -0.5 15572 2.25 0
15 373 7.2 -2.2 11509 1.62 0
16 404 7.5 -2.9 25447 1.00 0
17 385 7.7 -5.1 24090 0.66 0
18 365 7.8 -7.2 27786 0.31 0
19 366 7.7 -7.9 26195 0.25 0
20 421 7.7 -10.9 20516 0.25 0
21 354 7.8 -12.7 22759 0.25 0
22 367 8.0 -14.0 19028 0.25 0
23 413 8.1 -15.6 16971 0.25 0
24 362 8.1 -16.0 20036 0.25 0
25 385 8.0 -17.2 22485 0.25 1
26 343 8.1 -17.6 18730 0.25 0
27 369 8.2 -15.5 14538 0.25 0
28 363 8.4 -13.7 27561 0.25 0
29 318 8.5 -11.4 25985 0.25 0
30 393 8.5 -9.2 34670 0.25 0
31 325 8.5 -6.3 32066 0.25 0
32 403 8.5 -3.1 27186 0.25 0
33 392 8.5 0.0 29586 0.25 0
34 409 8.4 3.0 21359 0.25 0
35 485 8.3 5.4 21553 0.25 0
36 423 8.2 7.6 19573 0.25 0
37 428 8.1 9.7 24256 0.25 1
38 431 7.9 12.0 22380 0.25 0
39 416 7.6 11.6 16167 0.25 0
40 330 7.3 10.0 27297 0.25 0
41 314 7.1 10.8 28287 0.25 0
42 345 7.0 11.3 33474 0.39 0
43 365 7.1 10.1 28229 0.50 0
44 417 7.1 9.4 28785 0.50 0
45 356 7.1 9.6 25597 0.65 0
46 477 7.3 7.9 18130 0.75 0
47 423 7.3 7.3 20198 0.75 0
48 386 7.3 6.2 22849 0.75 0
49 390 7.2 4.9 23118 0.57 1
50 407 7.2 3.6 21925 0.36 0
51 398 7.1 2.9 20801 0.25 0
52 327 7.1 3.1 18785 0.25 0
53 304 7.1 1.7 20659 0.25 0
54 378 7.2 0.6 29367 0.25 0
55 311 7.3 -0.4 23992 0.25 0
56 376 7.4 -1.1 20645 0.25 0
57 340 7.4 -2.9 22356 0.08 0
58 383 7.5 -2.8 17902 0.00 0
59 467 7.4 -3.0 15879 0.00 0
60 439 7.4 -3.2 16963 0.00 0
M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 0 0 0 0 0 0 0 0 0 0
2 1 0 0 0 0 0 0 0 0 0
3 0 1 0 0 0 0 0 0 0 0
4 0 0 1 0 0 0 0 0 0 0
5 0 0 0 1 0 0 0 0 0 0
6 0 0 0 0 1 0 0 0 0 0
7 0 0 0 0 0 1 0 0 0 0
8 0 0 0 0 0 0 1 0 0 0
9 0 0 0 0 0 0 0 1 0 0
10 0 0 0 0 0 0 0 0 1 0
11 0 0 0 0 0 0 0 0 0 1
12 0 0 0 0 0 0 0 0 0 0
13 0 0 0 0 0 0 0 0 0 0
14 1 0 0 0 0 0 0 0 0 0
15 0 1 0 0 0 0 0 0 0 0
16 0 0 1 0 0 0 0 0 0 0
17 0 0 0 1 0 0 0 0 0 0
18 0 0 0 0 1 0 0 0 0 0
19 0 0 0 0 0 1 0 0 0 0
20 0 0 0 0 0 0 1 0 0 0
21 0 0 0 0 0 0 0 1 0 0
22 0 0 0 0 0 0 0 0 1 0
23 0 0 0 0 0 0 0 0 0 1
24 0 0 0 0 0 0 0 0 0 0
25 0 0 0 0 0 0 0 0 0 0
26 1 0 0 0 0 0 0 0 0 0
27 0 1 0 0 0 0 0 0 0 0
28 0 0 1 0 0 0 0 0 0 0
29 0 0 0 1 0 0 0 0 0 0
30 0 0 0 0 1 0 0 0 0 0
31 0 0 0 0 0 1 0 0 0 0
32 0 0 0 0 0 0 1 0 0 0
33 0 0 0 0 0 0 0 1 0 0
34 0 0 0 0 0 0 0 0 1 0
35 0 0 0 0 0 0 0 0 0 1
36 0 0 0 0 0 0 0 0 0 0
37 0 0 0 0 0 0 0 0 0 0
38 1 0 0 0 0 0 0 0 0 0
39 0 1 0 0 0 0 0 0 0 0
40 0 0 1 0 0 0 0 0 0 0
41 0 0 0 1 0 0 0 0 0 0
42 0 0 0 0 1 0 0 0 0 0
43 0 0 0 0 0 1 0 0 0 0
44 0 0 0 0 0 0 1 0 0 0
45 0 0 0 0 0 0 0 1 0 0
46 0 0 0 0 0 0 0 0 1 0
47 0 0 0 0 0 0 0 0 0 1
48 0 0 0 0 0 0 0 0 0 0
49 0 0 0 0 0 0 0 0 0 0
50 1 0 0 0 0 0 0 0 0 0
51 0 1 0 0 0 0 0 0 0 0
52 0 0 1 0 0 0 0 0 0 0
53 0 0 0 1 0 0 0 0 0 0
54 0 0 0 0 1 0 0 0 0 0
55 0 0 0 0 0 1 0 0 0 0
56 0 0 0 0 0 0 1 0 0 0
57 0 0 0 0 0 0 0 1 0 0
58 0 0 0 0 0 0 0 0 1 0
59 0 0 0 0 0 0 0 0 0 1
60 0 0 0 0 0 0 0 0 0 0
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) werkloosheidsgraad uitvoer personenwagens
237.971793 25.754226 1.885971 -0.002274
rentetarieven M1 M2 M3
8.295662 9.215727 13.535492 -1.574276
M4 M5 M6 M7
-31.133951 -46.608251 7.672382 -28.089551
M8 M9 M10 M11
19.211822 -4.034552 9.366294 49.191967
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-43.415 -21.742 -6.572 19.004 61.754
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 237.971793 92.729354 2.566 0.0138 *
werkloosheidsgraad 25.754226 13.195490 1.952 0.0574 .
uitvoer 1.885971 0.706613 2.669 0.0106 *
personenwagens -0.002274 0.001807 -1.258 0.2151
rentetarieven 8.295662 4.844641 1.712 0.0939 .
M1 9.215727 21.260260 0.433 0.6668
M2 13.535492 20.434122 0.662 0.5112
M3 -1.574276 22.040462 -0.071 0.9434
M4 -31.133951 23.545389 -1.322 0.1929
M5 -46.608251 23.457321 -1.987 0.0532 .
M6 7.672382 29.360571 0.261 0.7951
M7 -28.089551 26.150712 -1.074 0.2886
M8 19.211822 22.409938 0.857 0.3959
M9 -4.034552 22.907178 -0.176 0.8610
M10 9.366294 20.432556 0.458 0.6489
M11 49.191967 20.589727 2.389 0.0212 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 32.23 on 44 degrees of freedom
Multiple R-squared: 0.6011, Adjusted R-squared: 0.4651
F-statistic: 4.42 on 15 and 44 DF, p-value: 5.74e-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.7722579 0.4554842 0.22774211
[2,] 0.9109851 0.1780298 0.08901489
[3,] 0.8557358 0.2885284 0.14426419
[4,] 0.9224508 0.1550983 0.07754915
[5,] 0.9030140 0.1939720 0.09698599
[6,] 0.9109953 0.1780093 0.08900467
[7,] 0.8743374 0.2513251 0.12566256
[8,] 0.8779985 0.2440030 0.12200148
[9,] 0.8451636 0.3096727 0.15483637
[10,] 0.8049064 0.3901872 0.19509362
[11,] 0.7680468 0.4639063 0.23195317
[12,] 0.7965242 0.4069516 0.20347579
[13,] 0.7476399 0.5047202 0.25236010
[14,] 0.6517208 0.6965585 0.34827924
[15,] 0.6557353 0.6885294 0.34426470
[16,] 0.5606530 0.8786940 0.43934702
[17,] 0.5440328 0.9119345 0.45596723
[18,] 0.4280857 0.8561715 0.57191426
[19,] 0.3654811 0.7309623 0.63451885
[20,] 0.2726617 0.5453233 0.72733833
[21,] 0.1829038 0.3658076 0.81709618
[22,] 0.1855374 0.3710747 0.81446265
[23,] 0.2619900 0.5239800 0.73800999
> postscript(file="/var/fisher/rcomp/tmp/1opcm1356195914.ps",horizontal=F,onefile=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/fisher/rcomp/tmp/2uirr1356195914.ps",horizontal=F,onefile=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/fisher/rcomp/tmp/3d8c81356195914.ps",horizontal=F,onefile=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/fisher/rcomp/tmp/4ycw01356195914.ps",horizontal=F,onefile=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/fisher/rcomp/tmp/5fm1n1356195914.ps",horizontal=F,onefile=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 = 60
Frequency = 1
1 2 3 4 5 6 7
21.363782 -26.246754 35.152875 -33.660833 6.410749 3.684972 11.309447
8 9 10 11 12 13 14
-18.821490 57.838808 -17.612541 -8.752039 -40.110645 -37.579813 57.318499
15 16 17 18 19 20 21
-31.951996 59.033036 54.241003 -7.348153 30.189973 30.635268 -7.199558
22 23 24 25 26 27 28
-18.781941 -16.842086 -10.927433 13.263245 -43.414572 -18.371305 26.250655
29 30 31 32 33 34 35
-13.771246 22.544390 -21.083204 -7.514395 4.341887 -13.845588 20.818893
36 37 38 39 40 41 42
1.935604 6.981593 2.209842 -3.325031 -23.717406 -18.350265 -29.367175
43 44 45 46 47 48 49
13.245435 20.528311 -26.094797 61.753821 -26.238657 -5.945057 -4.028807
50 51 52 53 54 55 56
10.132985 18.495457 -27.905452 -28.530241 10.485966 -33.661651 -24.827694
57 58 59 60
-28.886340 -11.513752 31.013889 55.047531
> postscript(file="/var/fisher/rcomp/tmp/6sxzw1356195914.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> dum <- cbind(lag(myerror,k=1),myerror)
> dum
Time Series:
Start = 0
End = 60
Frequency = 1
lag(myerror, k = 1) myerror
0 21.363782 NA
1 -26.246754 21.363782
2 35.152875 -26.246754
3 -33.660833 35.152875
4 6.410749 -33.660833
5 3.684972 6.410749
6 11.309447 3.684972
7 -18.821490 11.309447
8 57.838808 -18.821490
9 -17.612541 57.838808
10 -8.752039 -17.612541
11 -40.110645 -8.752039
12 -37.579813 -40.110645
13 57.318499 -37.579813
14 -31.951996 57.318499
15 59.033036 -31.951996
16 54.241003 59.033036
17 -7.348153 54.241003
18 30.189973 -7.348153
19 30.635268 30.189973
20 -7.199558 30.635268
21 -18.781941 -7.199558
22 -16.842086 -18.781941
23 -10.927433 -16.842086
24 13.263245 -10.927433
25 -43.414572 13.263245
26 -18.371305 -43.414572
27 26.250655 -18.371305
28 -13.771246 26.250655
29 22.544390 -13.771246
30 -21.083204 22.544390
31 -7.514395 -21.083204
32 4.341887 -7.514395
33 -13.845588 4.341887
34 20.818893 -13.845588
35 1.935604 20.818893
36 6.981593 1.935604
37 2.209842 6.981593
38 -3.325031 2.209842
39 -23.717406 -3.325031
40 -18.350265 -23.717406
41 -29.367175 -18.350265
42 13.245435 -29.367175
43 20.528311 13.245435
44 -26.094797 20.528311
45 61.753821 -26.094797
46 -26.238657 61.753821
47 -5.945057 -26.238657
48 -4.028807 -5.945057
49 10.132985 -4.028807
50 18.495457 10.132985
51 -27.905452 18.495457
52 -28.530241 -27.905452
53 10.485966 -28.530241
54 -33.661651 10.485966
55 -24.827694 -33.661651
56 -28.886340 -24.827694
57 -11.513752 -28.886340
58 31.013889 -11.513752
59 55.047531 31.013889
60 NA 55.047531
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -26.246754 21.363782
[2,] 35.152875 -26.246754
[3,] -33.660833 35.152875
[4,] 6.410749 -33.660833
[5,] 3.684972 6.410749
[6,] 11.309447 3.684972
[7,] -18.821490 11.309447
[8,] 57.838808 -18.821490
[9,] -17.612541 57.838808
[10,] -8.752039 -17.612541
[11,] -40.110645 -8.752039
[12,] -37.579813 -40.110645
[13,] 57.318499 -37.579813
[14,] -31.951996 57.318499
[15,] 59.033036 -31.951996
[16,] 54.241003 59.033036
[17,] -7.348153 54.241003
[18,] 30.189973 -7.348153
[19,] 30.635268 30.189973
[20,] -7.199558 30.635268
[21,] -18.781941 -7.199558
[22,] -16.842086 -18.781941
[23,] -10.927433 -16.842086
[24,] 13.263245 -10.927433
[25,] -43.414572 13.263245
[26,] -18.371305 -43.414572
[27,] 26.250655 -18.371305
[28,] -13.771246 26.250655
[29,] 22.544390 -13.771246
[30,] -21.083204 22.544390
[31,] -7.514395 -21.083204
[32,] 4.341887 -7.514395
[33,] -13.845588 4.341887
[34,] 20.818893 -13.845588
[35,] 1.935604 20.818893
[36,] 6.981593 1.935604
[37,] 2.209842 6.981593
[38,] -3.325031 2.209842
[39,] -23.717406 -3.325031
[40,] -18.350265 -23.717406
[41,] -29.367175 -18.350265
[42,] 13.245435 -29.367175
[43,] 20.528311 13.245435
[44,] -26.094797 20.528311
[45,] 61.753821 -26.094797
[46,] -26.238657 61.753821
[47,] -5.945057 -26.238657
[48,] -4.028807 -5.945057
[49,] 10.132985 -4.028807
[50,] 18.495457 10.132985
[51,] -27.905452 18.495457
[52,] -28.530241 -27.905452
[53,] 10.485966 -28.530241
[54,] -33.661651 10.485966
[55,] -24.827694 -33.661651
[56,] -28.886340 -24.827694
[57,] -11.513752 -28.886340
[58,] 31.013889 -11.513752
[59,] 55.047531 31.013889
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -26.246754 21.363782
2 35.152875 -26.246754
3 -33.660833 35.152875
4 6.410749 -33.660833
5 3.684972 6.410749
6 11.309447 3.684972
7 -18.821490 11.309447
8 57.838808 -18.821490
9 -17.612541 57.838808
10 -8.752039 -17.612541
11 -40.110645 -8.752039
12 -37.579813 -40.110645
13 57.318499 -37.579813
14 -31.951996 57.318499
15 59.033036 -31.951996
16 54.241003 59.033036
17 -7.348153 54.241003
18 30.189973 -7.348153
19 30.635268 30.189973
20 -7.199558 30.635268
21 -18.781941 -7.199558
22 -16.842086 -18.781941
23 -10.927433 -16.842086
24 13.263245 -10.927433
25 -43.414572 13.263245
26 -18.371305 -43.414572
27 26.250655 -18.371305
28 -13.771246 26.250655
29 22.544390 -13.771246
30 -21.083204 22.544390
31 -7.514395 -21.083204
32 4.341887 -7.514395
33 -13.845588 4.341887
34 20.818893 -13.845588
35 1.935604 20.818893
36 6.981593 1.935604
37 2.209842 6.981593
38 -3.325031 2.209842
39 -23.717406 -3.325031
40 -18.350265 -23.717406
41 -29.367175 -18.350265
42 13.245435 -29.367175
43 20.528311 13.245435
44 -26.094797 20.528311
45 61.753821 -26.094797
46 -26.238657 61.753821
47 -5.945057 -26.238657
48 -4.028807 -5.945057
49 10.132985 -4.028807
50 18.495457 10.132985
51 -27.905452 18.495457
52 -28.530241 -27.905452
53 10.485966 -28.530241
54 -33.661651 10.485966
55 -24.827694 -33.661651
56 -28.886340 -24.827694
57 -11.513752 -28.886340
58 31.013889 -11.513752
59 55.047531 31.013889
> 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/fisher/rcomp/tmp/7vb661356195914.ps",horizontal=F,onefile=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/fisher/rcomp/tmp/8iv7n1356195914.ps",horizontal=F,onefile=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/fisher/rcomp/tmp/9akcs1356195914.ps",horizontal=F,onefile=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/fisher/rcomp/tmp/10fqvk1356195914.ps",horizontal=F,onefile=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/fisher/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/fisher/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/fisher/rcomp/tmp/116c661356195914.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/fisher/rcomp/tmp/12hrok1356195914.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/fisher/rcomp/tmp/13eyic1356195915.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/fisher/rcomp/tmp/14rl9k1356195915.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/fisher/rcomp/tmp/159c5t1356195915.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/fisher/rcomp/tmp/168k2f1356195915.tab")
+ }
>
> try(system("convert tmp/1opcm1356195914.ps tmp/1opcm1356195914.png",intern=TRUE))
character(0)
> try(system("convert tmp/2uirr1356195914.ps tmp/2uirr1356195914.png",intern=TRUE))
character(0)
> try(system("convert tmp/3d8c81356195914.ps tmp/3d8c81356195914.png",intern=TRUE))
character(0)
> try(system("convert tmp/4ycw01356195914.ps tmp/4ycw01356195914.png",intern=TRUE))
character(0)
> try(system("convert tmp/5fm1n1356195914.ps tmp/5fm1n1356195914.png",intern=TRUE))
character(0)
> try(system("convert tmp/6sxzw1356195914.ps tmp/6sxzw1356195914.png",intern=TRUE))
character(0)
> try(system("convert tmp/7vb661356195914.ps tmp/7vb661356195914.png",intern=TRUE))
character(0)
> try(system("convert tmp/8iv7n1356195914.ps tmp/8iv7n1356195914.png",intern=TRUE))
character(0)
> try(system("convert tmp/9akcs1356195914.ps tmp/9akcs1356195914.png",intern=TRUE))
character(0)
> try(system("convert tmp/10fqvk1356195914.ps tmp/10fqvk1356195914.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
6.405 1.830 8.231