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(180144
+ ,11554.5
+ ,173666
+ ,13182.1
+ ,165688
+ ,14800.1
+ ,161570
+ ,12150.7
+ ,156145
+ ,14478.2
+ ,153730
+ ,13253.9
+ ,182698
+ ,12036.8
+ ,200765
+ ,12653.2
+ ,176512
+ ,14035.4
+ ,166618
+ ,14571.4
+ ,158644
+ ,15400.9
+ ,159585
+ ,14283.2
+ ,163095
+ ,14485.3
+ ,159044
+ ,14196.3
+ ,155511
+ ,15559.1
+ ,153745
+ ,13767.4
+ ,150569
+ ,14634
+ ,150605
+ ,14381.1
+ ,179612
+ ,12509.9
+ ,194690
+ ,12122.3
+ ,189917
+ ,13122.3
+ ,184128
+ ,13908.7
+ ,175335
+ ,13456.5
+ ,179566
+ ,12441.6
+ ,181140
+ ,12953
+ ,177876
+ ,13057.2
+ ,175041
+ ,14350.1
+ ,169292
+ ,13830.2
+ ,166070
+ ,13755.5
+ ,166972
+ ,13574.4
+ ,206348
+ ,12802.6
+ ,215706
+ ,11737.3
+ ,202108
+ ,13850.2
+ ,195411
+ ,15081.8
+ ,193111
+ ,13653.3
+ ,195198
+ ,14019.1
+ ,198770
+ ,13962
+ ,194163
+ ,13768.7
+ ,190420
+ ,14747.1
+ ,189733
+ ,13858.1
+ ,186029
+ ,13188
+ ,191531
+ ,13693.1
+ ,232571
+ ,12970
+ ,243477
+ ,11392.8
+ ,227247
+ ,13985.2
+ ,217859
+ ,14994.7
+ ,208679
+ ,13584.7
+ ,213188
+ ,14257.8
+ ,216234
+ ,13553.4
+ ,213586
+ ,14007.3
+ ,209465
+ ,16535.8
+ ,204045
+ ,14721.4
+ ,200237
+ ,13664.6
+ ,203666
+ ,16405.9
+ ,241476
+ ,13829.4
+ ,260307
+ ,13735.6
+ ,243324
+ ,15870.5
+ ,244460
+ ,15962.4
+ ,233575
+ ,15744.1
+ ,237217
+ ,16083.7
+ ,235243
+ ,14863.9
+ ,230354
+ ,15533.1
+ ,227184
+ ,17473.1
+ ,221678
+ ,15925.5
+ ,217142
+ ,15573.7
+ ,219452
+ ,17495
+ ,256446
+ ,14155.8
+ ,265845
+ ,14913.9
+ ,248624
+ ,17250.4
+ ,241114
+ ,15879.8
+ ,229245
+ ,17647.8
+ ,231805
+ ,17749.9
+ ,219277
+ ,17111.8
+ ,219313
+ ,16934.8
+ ,212610
+ ,20280
+ ,214771
+ ,16238.2
+ ,211142
+ ,17896.1
+ ,211457
+ ,18089.3
+ ,240048
+ ,15660
+ ,240636
+ ,16162.4
+ ,230580
+ ,17850.1
+ ,208795
+ ,18520.4
+ ,197922
+ ,18524.7
+ ,194596
+ ,16843.7)
+ ,dim=c(2
+ ,84)
+ ,dimnames=list(c('werkloosheid'
+ ,'invoer')
+ ,1:84))
> y <- array(NA,dim=c(2,84),dimnames=list(c('werkloosheid','invoer'),1:84))
> 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
werkloosheid invoer M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 180144 11554.5 1 0 0 0 0 0 0 0 0 0 0 1
2 173666 13182.1 0 1 0 0 0 0 0 0 0 0 0 2
3 165688 14800.1 0 0 1 0 0 0 0 0 0 0 0 3
4 161570 12150.7 0 0 0 1 0 0 0 0 0 0 0 4
5 156145 14478.2 0 0 0 0 1 0 0 0 0 0 0 5
6 153730 13253.9 0 0 0 0 0 1 0 0 0 0 0 6
7 182698 12036.8 0 0 0 0 0 0 1 0 0 0 0 7
8 200765 12653.2 0 0 0 0 0 0 0 1 0 0 0 8
9 176512 14035.4 0 0 0 0 0 0 0 0 1 0 0 9
10 166618 14571.4 0 0 0 0 0 0 0 0 0 1 0 10
11 158644 15400.9 0 0 0 0 0 0 0 0 0 0 1 11
12 159585 14283.2 0 0 0 0 0 0 0 0 0 0 0 12
13 163095 14485.3 1 0 0 0 0 0 0 0 0 0 0 13
14 159044 14196.3 0 1 0 0 0 0 0 0 0 0 0 14
15 155511 15559.1 0 0 1 0 0 0 0 0 0 0 0 15
16 153745 13767.4 0 0 0 1 0 0 0 0 0 0 0 16
17 150569 14634.0 0 0 0 0 1 0 0 0 0 0 0 17
18 150605 14381.1 0 0 0 0 0 1 0 0 0 0 0 18
19 179612 12509.9 0 0 0 0 0 0 1 0 0 0 0 19
20 194690 12122.3 0 0 0 0 0 0 0 1 0 0 0 20
21 189917 13122.3 0 0 0 0 0 0 0 0 1 0 0 21
22 184128 13908.7 0 0 0 0 0 0 0 0 0 1 0 22
23 175335 13456.5 0 0 0 0 0 0 0 0 0 0 1 23
24 179566 12441.6 0 0 0 0 0 0 0 0 0 0 0 24
25 181140 12953.0 1 0 0 0 0 0 0 0 0 0 0 25
26 177876 13057.2 0 1 0 0 0 0 0 0 0 0 0 26
27 175041 14350.1 0 0 1 0 0 0 0 0 0 0 0 27
28 169292 13830.2 0 0 0 1 0 0 0 0 0 0 0 28
29 166070 13755.5 0 0 0 0 1 0 0 0 0 0 0 29
30 166972 13574.4 0 0 0 0 0 1 0 0 0 0 0 30
31 206348 12802.6 0 0 0 0 0 0 1 0 0 0 0 31
32 215706 11737.3 0 0 0 0 0 0 0 1 0 0 0 32
33 202108 13850.2 0 0 0 0 0 0 0 0 1 0 0 33
34 195411 15081.8 0 0 0 0 0 0 0 0 0 1 0 34
35 193111 13653.3 0 0 0 0 0 0 0 0 0 0 1 35
36 195198 14019.1 0 0 0 0 0 0 0 0 0 0 0 36
37 198770 13962.0 1 0 0 0 0 0 0 0 0 0 0 37
38 194163 13768.7 0 1 0 0 0 0 0 0 0 0 0 38
39 190420 14747.1 0 0 1 0 0 0 0 0 0 0 0 39
40 189733 13858.1 0 0 0 1 0 0 0 0 0 0 0 40
41 186029 13188.0 0 0 0 0 1 0 0 0 0 0 0 41
42 191531 13693.1 0 0 0 0 0 1 0 0 0 0 0 42
43 232571 12970.0 0 0 0 0 0 0 1 0 0 0 0 43
44 243477 11392.8 0 0 0 0 0 0 0 1 0 0 0 44
45 227247 13985.2 0 0 0 0 0 0 0 0 1 0 0 45
46 217859 14994.7 0 0 0 0 0 0 0 0 0 1 0 46
47 208679 13584.7 0 0 0 0 0 0 0 0 0 0 1 47
48 213188 14257.8 0 0 0 0 0 0 0 0 0 0 0 48
49 216234 13553.4 1 0 0 0 0 0 0 0 0 0 0 49
50 213586 14007.3 0 1 0 0 0 0 0 0 0 0 0 50
51 209465 16535.8 0 0 1 0 0 0 0 0 0 0 0 51
52 204045 14721.4 0 0 0 1 0 0 0 0 0 0 0 52
53 200237 13664.6 0 0 0 0 1 0 0 0 0 0 0 53
54 203666 16405.9 0 0 0 0 0 1 0 0 0 0 0 54
55 241476 13829.4 0 0 0 0 0 0 1 0 0 0 0 55
56 260307 13735.6 0 0 0 0 0 0 0 1 0 0 0 56
57 243324 15870.5 0 0 0 0 0 0 0 0 1 0 0 57
58 244460 15962.4 0 0 0 0 0 0 0 0 0 1 0 58
59 233575 15744.1 0 0 0 0 0 0 0 0 0 0 1 59
60 237217 16083.7 0 0 0 0 0 0 0 0 0 0 0 60
61 235243 14863.9 1 0 0 0 0 0 0 0 0 0 0 61
62 230354 15533.1 0 1 0 0 0 0 0 0 0 0 0 62
63 227184 17473.1 0 0 1 0 0 0 0 0 0 0 0 63
64 221678 15925.5 0 0 0 1 0 0 0 0 0 0 0 64
65 217142 15573.7 0 0 0 0 1 0 0 0 0 0 0 65
66 219452 17495.0 0 0 0 0 0 1 0 0 0 0 0 66
67 256446 14155.8 0 0 0 0 0 0 1 0 0 0 0 67
68 265845 14913.9 0 0 0 0 0 0 0 1 0 0 0 68
69 248624 17250.4 0 0 0 0 0 0 0 0 1 0 0 69
70 241114 15879.8 0 0 0 0 0 0 0 0 0 1 0 70
71 229245 17647.8 0 0 0 0 0 0 0 0 0 0 1 71
72 231805 17749.9 0 0 0 0 0 0 0 0 0 0 0 72
73 219277 17111.8 1 0 0 0 0 0 0 0 0 0 0 73
74 219313 16934.8 0 1 0 0 0 0 0 0 0 0 0 74
75 212610 20280.0 0 0 1 0 0 0 0 0 0 0 0 75
76 214771 16238.2 0 0 0 1 0 0 0 0 0 0 0 76
77 211142 17896.1 0 0 0 0 1 0 0 0 0 0 0 77
78 211457 18089.3 0 0 0 0 0 1 0 0 0 0 0 78
79 240048 15660.0 0 0 0 0 0 0 1 0 0 0 0 79
80 240636 16162.4 0 0 0 0 0 0 0 1 0 0 0 80
81 230580 17850.1 0 0 0 0 0 0 0 0 1 0 0 81
82 208795 18520.4 0 0 0 0 0 0 0 0 0 1 0 82
83 197922 18524.7 0 0 0 0 0 0 0 0 0 0 1 83
84 194596 16843.7 0 0 0 0 0 0 0 0 0 0 0 84
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) invoer M1 M2 M3 M4
220056.614 -5.046 5574.863 2255.008 5887.486 -7880.629
M5 M6 M7 M8 M9 M10
-11066.338 -8159.914 15859.671 25505.174 19120.690 11487.443
M11 t
792.024 1202.347
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-41469 -7093 -907 9478 26174
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 220056.614 19897.546 11.059 < 2e-16 ***
invoer -5.046 1.531 -3.297 0.001539 **
M1 5574.863 6921.634 0.805 0.423302
M2 2255.008 6890.339 0.327 0.744440
M3 5887.486 7326.610 0.804 0.424362
M4 -7880.629 6892.360 -1.143 0.256774
M5 -11066.338 6873.056 -1.610 0.111876
M6 -8159.914 6912.424 -1.180 0.241809
M7 15859.671 7193.680 2.205 0.030766 *
M8 25505.174 7305.661 3.491 0.000837 ***
M9 19120.690 6870.163 2.783 0.006916 **
M10 11487.443 6917.440 1.661 0.101257
M11 792.024 6886.678 0.115 0.908768
t 1202.347 101.986 11.789 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 12840 on 70 degrees of freedom
Multiple R-squared: 0.8398, Adjusted R-squared: 0.8101
F-statistic: 28.23 on 13 and 70 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,] 9.707768e-03 1.941554e-02 0.9902922
[2,] 6.134055e-03 1.226811e-02 0.9938659
[3,] 1.724896e-03 3.449791e-03 0.9982751
[4,] 4.581129e-04 9.162258e-04 0.9995419
[5,] 1.751368e-03 3.502736e-03 0.9982486
[6,] 3.273066e-03 6.546131e-03 0.9967269
[7,] 1.136737e-03 2.273475e-03 0.9988633
[8,] 4.523210e-04 9.046419e-04 0.9995477
[9,] 2.384480e-04 4.768961e-04 0.9997616
[10,] 8.684552e-05 1.736910e-04 0.9999132
[11,] 3.435077e-05 6.870154e-05 0.9999656
[12,] 9.644124e-05 1.928825e-04 0.9999036
[13,] 3.651476e-05 7.302953e-05 0.9999635
[14,] 2.559858e-05 5.119715e-05 0.9999744
[15,] 4.301258e-04 8.602516e-04 0.9995699
[16,] 2.706043e-04 5.412086e-04 0.9997294
[17,] 3.377167e-04 6.754334e-04 0.9996623
[18,] 6.881516e-04 1.376303e-03 0.9993118
[19,] 5.632708e-04 1.126542e-03 0.9994367
[20,] 1.229111e-03 2.458222e-03 0.9987709
[21,] 1.563378e-03 3.126755e-03 0.9984366
[22,] 1.573991e-03 3.147982e-03 0.9984260
[23,] 1.232373e-03 2.464745e-03 0.9987676
[24,] 1.875186e-03 3.750372e-03 0.9981248
[25,] 1.700131e-03 3.400261e-03 0.9982999
[26,] 2.220283e-03 4.440565e-03 0.9977797
[27,] 9.415526e-03 1.883105e-02 0.9905845
[28,] 8.722133e-03 1.744427e-02 0.9912779
[29,] 1.123767e-02 2.247535e-02 0.9887623
[30,] 1.664402e-02 3.328805e-02 0.9833560
[31,] 1.298127e-02 2.596255e-02 0.9870187
[32,] 1.476057e-02 2.952115e-02 0.9852394
[33,] 1.293424e-02 2.586848e-02 0.9870658
[34,] 1.342628e-02 2.685255e-02 0.9865737
[35,] 1.779216e-02 3.558432e-02 0.9822078
[36,] 3.391435e-02 6.782869e-02 0.9660857
[37,] 6.808337e-02 1.361667e-01 0.9319166
[38,] 2.010434e-01 4.020869e-01 0.7989566
[39,] 4.580145e-01 9.160291e-01 0.5419855
[40,] 5.528726e-01 8.942548e-01 0.4471274
[41,] 7.671875e-01 4.656249e-01 0.2328125
[42,] 7.791777e-01 4.416446e-01 0.2208223
[43,] 7.491441e-01 5.017117e-01 0.2508559
[44,] 7.082320e-01 5.835360e-01 0.2917680
[45,] 6.114849e-01 7.770302e-01 0.3885151
[46,] 5.403826e-01 9.192347e-01 0.4596174
[47,] 4.281503e-01 8.563006e-01 0.5718497
[48,] 4.821034e-01 9.642067e-01 0.5178966
[49,] 5.066947e-01 9.866106e-01 0.4933053
[50,] 6.966467e-01 6.067067e-01 0.3033533
[51,] 6.877570e-01 6.244860e-01 0.3122430
> postscript(file="/var/www/html/freestat/rcomp/tmp/1u0ji1229718605.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/freestat/rcomp/tmp/272tc1229718605.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/freestat/rcomp/tmp/3kklk1229718605.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/freestat/rcomp/tmp/4if2a1229718605.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/freestat/rcomp/tmp/5kevz1229718605.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 = 84
Frequency = 1
1 2 3 4 5 6
11611.3030 15463.2740 10814.4744 5894.0297 14196.3770 1495.0934
7 8 9 10 11 12
-900.0220 9429.3289 -2667.2975 -3425.8748 2278.6475 -2830.3114
13 14 15 16 17 18
-5077.7747 -8469.4896 -9960.9648 -8201.6691 -5021.6590 -10370.5003
19 20 21 22 23 24
-16027.0413 -13752.6241 -8297.7365 -3687.8578 -5269.4741 -6569.7299
25 26 27 28 29 30
-9192.5426 -9813.2682 -10959.4414 -6765.9599 -8381.5150 -12502.0714
31 32 33 34 35 36
-2242.3141 -9107.4021 -6862.0985 -913.8516 -928.6344 2593.7779
37 38 39 40 41 42
-899.5441 -4364.3806 -8005.4423 -612.3474 -5714.1426 -1772.3048
43 44 45 46 47 48
10397.1806 2497.1727 4529.9138 6666.4995 -134.9369 7360.0346
49 50 51 52 53 54
74.5982 1834.3716 5636.7280 3627.4852 -3529.5018 9622.6438
55 56 57 58 59 60
9210.3348 16720.1935 15691.5037 23722.1085 21228.6932 26173.9068
61 62 63 64 65 66
11267.8905 11873.0140 13656.9462 12907.9096 8580.1770 16475.8069
67 68 69 70 71 72
11399.1038 13775.4376 13525.9711 5531.1655 12076.1249 14740.9727
73 74 75 76 77 78
-7783.9304 -6523.5212 -1182.3001 -6849.4481 -129.7355 -2948.6675
79 80 81 82 83 84
-11837.2419 -19562.1066 -15920.2562 -27892.1892 -29250.4203 -41468.6506
> postscript(file="/var/www/html/freestat/rcomp/tmp/6wiuk1229718605.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 = 84
Frequency = 1
lag(myerror, k = 1) myerror
0 11611.3030 NA
1 15463.2740 11611.3030
2 10814.4744 15463.2740
3 5894.0297 10814.4744
4 14196.3770 5894.0297
5 1495.0934 14196.3770
6 -900.0220 1495.0934
7 9429.3289 -900.0220
8 -2667.2975 9429.3289
9 -3425.8748 -2667.2975
10 2278.6475 -3425.8748
11 -2830.3114 2278.6475
12 -5077.7747 -2830.3114
13 -8469.4896 -5077.7747
14 -9960.9648 -8469.4896
15 -8201.6691 -9960.9648
16 -5021.6590 -8201.6691
17 -10370.5003 -5021.6590
18 -16027.0413 -10370.5003
19 -13752.6241 -16027.0413
20 -8297.7365 -13752.6241
21 -3687.8578 -8297.7365
22 -5269.4741 -3687.8578
23 -6569.7299 -5269.4741
24 -9192.5426 -6569.7299
25 -9813.2682 -9192.5426
26 -10959.4414 -9813.2682
27 -6765.9599 -10959.4414
28 -8381.5150 -6765.9599
29 -12502.0714 -8381.5150
30 -2242.3141 -12502.0714
31 -9107.4021 -2242.3141
32 -6862.0985 -9107.4021
33 -913.8516 -6862.0985
34 -928.6344 -913.8516
35 2593.7779 -928.6344
36 -899.5441 2593.7779
37 -4364.3806 -899.5441
38 -8005.4423 -4364.3806
39 -612.3474 -8005.4423
40 -5714.1426 -612.3474
41 -1772.3048 -5714.1426
42 10397.1806 -1772.3048
43 2497.1727 10397.1806
44 4529.9138 2497.1727
45 6666.4995 4529.9138
46 -134.9369 6666.4995
47 7360.0346 -134.9369
48 74.5982 7360.0346
49 1834.3716 74.5982
50 5636.7280 1834.3716
51 3627.4852 5636.7280
52 -3529.5018 3627.4852
53 9622.6438 -3529.5018
54 9210.3348 9622.6438
55 16720.1935 9210.3348
56 15691.5037 16720.1935
57 23722.1085 15691.5037
58 21228.6932 23722.1085
59 26173.9068 21228.6932
60 11267.8905 26173.9068
61 11873.0140 11267.8905
62 13656.9462 11873.0140
63 12907.9096 13656.9462
64 8580.1770 12907.9096
65 16475.8069 8580.1770
66 11399.1038 16475.8069
67 13775.4376 11399.1038
68 13525.9711 13775.4376
69 5531.1655 13525.9711
70 12076.1249 5531.1655
71 14740.9727 12076.1249
72 -7783.9304 14740.9727
73 -6523.5212 -7783.9304
74 -1182.3001 -6523.5212
75 -6849.4481 -1182.3001
76 -129.7355 -6849.4481
77 -2948.6675 -129.7355
78 -11837.2419 -2948.6675
79 -19562.1066 -11837.2419
80 -15920.2562 -19562.1066
81 -27892.1892 -15920.2562
82 -29250.4203 -27892.1892
83 -41468.6506 -29250.4203
84 NA -41468.6506
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 15463.2740 11611.3030
[2,] 10814.4744 15463.2740
[3,] 5894.0297 10814.4744
[4,] 14196.3770 5894.0297
[5,] 1495.0934 14196.3770
[6,] -900.0220 1495.0934
[7,] 9429.3289 -900.0220
[8,] -2667.2975 9429.3289
[9,] -3425.8748 -2667.2975
[10,] 2278.6475 -3425.8748
[11,] -2830.3114 2278.6475
[12,] -5077.7747 -2830.3114
[13,] -8469.4896 -5077.7747
[14,] -9960.9648 -8469.4896
[15,] -8201.6691 -9960.9648
[16,] -5021.6590 -8201.6691
[17,] -10370.5003 -5021.6590
[18,] -16027.0413 -10370.5003
[19,] -13752.6241 -16027.0413
[20,] -8297.7365 -13752.6241
[21,] -3687.8578 -8297.7365
[22,] -5269.4741 -3687.8578
[23,] -6569.7299 -5269.4741
[24,] -9192.5426 -6569.7299
[25,] -9813.2682 -9192.5426
[26,] -10959.4414 -9813.2682
[27,] -6765.9599 -10959.4414
[28,] -8381.5150 -6765.9599
[29,] -12502.0714 -8381.5150
[30,] -2242.3141 -12502.0714
[31,] -9107.4021 -2242.3141
[32,] -6862.0985 -9107.4021
[33,] -913.8516 -6862.0985
[34,] -928.6344 -913.8516
[35,] 2593.7779 -928.6344
[36,] -899.5441 2593.7779
[37,] -4364.3806 -899.5441
[38,] -8005.4423 -4364.3806
[39,] -612.3474 -8005.4423
[40,] -5714.1426 -612.3474
[41,] -1772.3048 -5714.1426
[42,] 10397.1806 -1772.3048
[43,] 2497.1727 10397.1806
[44,] 4529.9138 2497.1727
[45,] 6666.4995 4529.9138
[46,] -134.9369 6666.4995
[47,] 7360.0346 -134.9369
[48,] 74.5982 7360.0346
[49,] 1834.3716 74.5982
[50,] 5636.7280 1834.3716
[51,] 3627.4852 5636.7280
[52,] -3529.5018 3627.4852
[53,] 9622.6438 -3529.5018
[54,] 9210.3348 9622.6438
[55,] 16720.1935 9210.3348
[56,] 15691.5037 16720.1935
[57,] 23722.1085 15691.5037
[58,] 21228.6932 23722.1085
[59,] 26173.9068 21228.6932
[60,] 11267.8905 26173.9068
[61,] 11873.0140 11267.8905
[62,] 13656.9462 11873.0140
[63,] 12907.9096 13656.9462
[64,] 8580.1770 12907.9096
[65,] 16475.8069 8580.1770
[66,] 11399.1038 16475.8069
[67,] 13775.4376 11399.1038
[68,] 13525.9711 13775.4376
[69,] 5531.1655 13525.9711
[70,] 12076.1249 5531.1655
[71,] 14740.9727 12076.1249
[72,] -7783.9304 14740.9727
[73,] -6523.5212 -7783.9304
[74,] -1182.3001 -6523.5212
[75,] -6849.4481 -1182.3001
[76,] -129.7355 -6849.4481
[77,] -2948.6675 -129.7355
[78,] -11837.2419 -2948.6675
[79,] -19562.1066 -11837.2419
[80,] -15920.2562 -19562.1066
[81,] -27892.1892 -15920.2562
[82,] -29250.4203 -27892.1892
[83,] -41468.6506 -29250.4203
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 15463.2740 11611.3030
2 10814.4744 15463.2740
3 5894.0297 10814.4744
4 14196.3770 5894.0297
5 1495.0934 14196.3770
6 -900.0220 1495.0934
7 9429.3289 -900.0220
8 -2667.2975 9429.3289
9 -3425.8748 -2667.2975
10 2278.6475 -3425.8748
11 -2830.3114 2278.6475
12 -5077.7747 -2830.3114
13 -8469.4896 -5077.7747
14 -9960.9648 -8469.4896
15 -8201.6691 -9960.9648
16 -5021.6590 -8201.6691
17 -10370.5003 -5021.6590
18 -16027.0413 -10370.5003
19 -13752.6241 -16027.0413
20 -8297.7365 -13752.6241
21 -3687.8578 -8297.7365
22 -5269.4741 -3687.8578
23 -6569.7299 -5269.4741
24 -9192.5426 -6569.7299
25 -9813.2682 -9192.5426
26 -10959.4414 -9813.2682
27 -6765.9599 -10959.4414
28 -8381.5150 -6765.9599
29 -12502.0714 -8381.5150
30 -2242.3141 -12502.0714
31 -9107.4021 -2242.3141
32 -6862.0985 -9107.4021
33 -913.8516 -6862.0985
34 -928.6344 -913.8516
35 2593.7779 -928.6344
36 -899.5441 2593.7779
37 -4364.3806 -899.5441
38 -8005.4423 -4364.3806
39 -612.3474 -8005.4423
40 -5714.1426 -612.3474
41 -1772.3048 -5714.1426
42 10397.1806 -1772.3048
43 2497.1727 10397.1806
44 4529.9138 2497.1727
45 6666.4995 4529.9138
46 -134.9369 6666.4995
47 7360.0346 -134.9369
48 74.5982 7360.0346
49 1834.3716 74.5982
50 5636.7280 1834.3716
51 3627.4852 5636.7280
52 -3529.5018 3627.4852
53 9622.6438 -3529.5018
54 9210.3348 9622.6438
55 16720.1935 9210.3348
56 15691.5037 16720.1935
57 23722.1085 15691.5037
58 21228.6932 23722.1085
59 26173.9068 21228.6932
60 11267.8905 26173.9068
61 11873.0140 11267.8905
62 13656.9462 11873.0140
63 12907.9096 13656.9462
64 8580.1770 12907.9096
65 16475.8069 8580.1770
66 11399.1038 16475.8069
67 13775.4376 11399.1038
68 13525.9711 13775.4376
69 5531.1655 13525.9711
70 12076.1249 5531.1655
71 14740.9727 12076.1249
72 -7783.9304 14740.9727
73 -6523.5212 -7783.9304
74 -1182.3001 -6523.5212
75 -6849.4481 -1182.3001
76 -129.7355 -6849.4481
77 -2948.6675 -129.7355
78 -11837.2419 -2948.6675
79 -19562.1066 -11837.2419
80 -15920.2562 -19562.1066
81 -27892.1892 -15920.2562
82 -29250.4203 -27892.1892
83 -41468.6506 -29250.4203
> 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/freestat/rcomp/tmp/7l3wd1229718605.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/freestat/rcomp/tmp/8fknl1229718605.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/freestat/rcomp/tmp/9cq0t1229718605.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/freestat/rcomp/tmp/10fb1y1229718605.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/freestat/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/html/freestat/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/freestat/rcomp/tmp/11it6n1229718605.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/freestat/rcomp/tmp/12kb741229718606.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/freestat/rcomp/tmp/13f8fs1229718606.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/freestat/rcomp/tmp/14a7qn1229718606.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/freestat/rcomp/tmp/156hjc1229718606.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/freestat/rcomp/tmp/166sl91229718606.tab")
+ }
>
> system("convert tmp/1u0ji1229718605.ps tmp/1u0ji1229718605.png")
> system("convert tmp/272tc1229718605.ps tmp/272tc1229718605.png")
> system("convert tmp/3kklk1229718605.ps tmp/3kklk1229718605.png")
> system("convert tmp/4if2a1229718605.ps tmp/4if2a1229718605.png")
> system("convert tmp/5kevz1229718605.ps tmp/5kevz1229718605.png")
> system("convert tmp/6wiuk1229718605.ps tmp/6wiuk1229718605.png")
> system("convert tmp/7l3wd1229718605.ps tmp/7l3wd1229718605.png")
> system("convert tmp/8fknl1229718605.ps tmp/8fknl1229718605.png")
> system("convert tmp/9cq0t1229718605.ps tmp/9cq0t1229718605.png")
> system("convert tmp/10fb1y1229718605.ps tmp/10fb1y1229718605.png")
>
>
> proc.time()
user system elapsed
4.034 2.551 4.369