R version 2.7.0 (2008-04-22)
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(-3.3,0,-3.5,0,-3.5,0,-8.4,0,-15.7,0,-18.7,0,-22.8,0,-20.7,0,-14,0,-6.3,0,0.7,1,0.2,1,0.8,1,1.2,1,4.5,1,0.4,1,5.9,1,6.5,1,12.8,1,4.2,1,-3.3,0,-12.5,0,-16.3,0,-10.5,0,-11.8,0,-11.4,0,-17.7,0,-17.3,0,-18.6,0,-17.9,0,-21.4,0,-19.4,0,-15.5,0,-7.7,0,-0.7,0,-1.6,0,1.4,1,0.7,1,9.5,1,1.4,1,4.1,1,6.6,1,18.4,1,16.9,1,9.2,1,-4.3,0,-5.9,0,-7.7,0,-5.4,0,-2.3,0,-4.8,0,2.3,0,-5.2,0,-10,0,-17.1,0,-14.4,0,-3.9,0,3.7,1,6.5,1,0.9,1,-4.1,0,-7,0,-12.2,0,-2.5,0,4.4,1,13.7,1,12.3,1,13.4,1,2.2,1,1.7,1,-7.2,0,-4.8,0,-2.9,0,-2.4,0,-2.5,0,-5.3,0,-7.1,0,-8,0,-8.9,0,-7.7,0,-1.1,0,4,1,9.6,1,10.9,1,13,1,14.9,1,20.1,1,10.8,1,11,1,3.8,1,10.8,1,7.6,1,10.2,1,2.2,1,-0.1,0,-1.7,0,-4.8,0),dim=c(2,97),dimnames=list(c('Registraties','D'),1:97))
> y <- array(NA,dim=c(2,97),dimnames=list(c('Registraties','D'),1:97))
> 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
Registraties D M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 -3.3 0 1 0 0 0 0 0 0 0 0 0 0 1
2 -3.5 0 0 1 0 0 0 0 0 0 0 0 0 2
3 -3.5 0 0 0 1 0 0 0 0 0 0 0 0 3
4 -8.4 0 0 0 0 1 0 0 0 0 0 0 0 4
5 -15.7 0 0 0 0 0 1 0 0 0 0 0 0 5
6 -18.7 0 0 0 0 0 0 1 0 0 0 0 0 6
7 -22.8 0 0 0 0 0 0 0 1 0 0 0 0 7
8 -20.7 0 0 0 0 0 0 0 0 1 0 0 0 8
9 -14.0 0 0 0 0 0 0 0 0 0 1 0 0 9
10 -6.3 0 0 0 0 0 0 0 0 0 0 1 0 10
11 0.7 1 0 0 0 0 0 0 0 0 0 0 1 11
12 0.2 1 0 0 0 0 0 0 0 0 0 0 0 12
13 0.8 1 1 0 0 0 0 0 0 0 0 0 0 13
14 1.2 1 0 1 0 0 0 0 0 0 0 0 0 14
15 4.5 1 0 0 1 0 0 0 0 0 0 0 0 15
16 0.4 1 0 0 0 1 0 0 0 0 0 0 0 16
17 5.9 1 0 0 0 0 1 0 0 0 0 0 0 17
18 6.5 1 0 0 0 0 0 1 0 0 0 0 0 18
19 12.8 1 0 0 0 0 0 0 1 0 0 0 0 19
20 4.2 1 0 0 0 0 0 0 0 1 0 0 0 20
21 -3.3 0 0 0 0 0 0 0 0 0 1 0 0 21
22 -12.5 0 0 0 0 0 0 0 0 0 0 1 0 22
23 -16.3 0 0 0 0 0 0 0 0 0 0 0 1 23
24 -10.5 0 0 0 0 0 0 0 0 0 0 0 0 24
25 -11.8 0 1 0 0 0 0 0 0 0 0 0 0 25
26 -11.4 0 0 1 0 0 0 0 0 0 0 0 0 26
27 -17.7 0 0 0 1 0 0 0 0 0 0 0 0 27
28 -17.3 0 0 0 0 1 0 0 0 0 0 0 0 28
29 -18.6 0 0 0 0 0 1 0 0 0 0 0 0 29
30 -17.9 0 0 0 0 0 0 1 0 0 0 0 0 30
31 -21.4 0 0 0 0 0 0 0 1 0 0 0 0 31
32 -19.4 0 0 0 0 0 0 0 0 1 0 0 0 32
33 -15.5 0 0 0 0 0 0 0 0 0 1 0 0 33
34 -7.7 0 0 0 0 0 0 0 0 0 0 1 0 34
35 -0.7 0 0 0 0 0 0 0 0 0 0 0 1 35
36 -1.6 0 0 0 0 0 0 0 0 0 0 0 0 36
37 1.4 1 1 0 0 0 0 0 0 0 0 0 0 37
38 0.7 1 0 1 0 0 0 0 0 0 0 0 0 38
39 9.5 1 0 0 1 0 0 0 0 0 0 0 0 39
40 1.4 1 0 0 0 1 0 0 0 0 0 0 0 40
41 4.1 1 0 0 0 0 1 0 0 0 0 0 0 41
42 6.6 1 0 0 0 0 0 1 0 0 0 0 0 42
43 18.4 1 0 0 0 0 0 0 1 0 0 0 0 43
44 16.9 1 0 0 0 0 0 0 0 1 0 0 0 44
45 9.2 1 0 0 0 0 0 0 0 0 1 0 0 45
46 -4.3 0 0 0 0 0 0 0 0 0 0 1 0 46
47 -5.9 0 0 0 0 0 0 0 0 0 0 0 1 47
48 -7.7 0 0 0 0 0 0 0 0 0 0 0 0 48
49 -5.4 0 1 0 0 0 0 0 0 0 0 0 0 49
50 -2.3 0 0 1 0 0 0 0 0 0 0 0 0 50
51 -4.8 0 0 0 1 0 0 0 0 0 0 0 0 51
52 2.3 0 0 0 0 1 0 0 0 0 0 0 0 52
53 -5.2 0 0 0 0 0 1 0 0 0 0 0 0 53
54 -10.0 0 0 0 0 0 0 1 0 0 0 0 0 54
55 -17.1 0 0 0 0 0 0 0 1 0 0 0 0 55
56 -14.4 0 0 0 0 0 0 0 0 1 0 0 0 56
57 -3.9 0 0 0 0 0 0 0 0 0 1 0 0 57
58 3.7 1 0 0 0 0 0 0 0 0 0 1 0 58
59 6.5 1 0 0 0 0 0 0 0 0 0 0 1 59
60 0.9 1 0 0 0 0 0 0 0 0 0 0 0 60
61 -4.1 0 1 0 0 0 0 0 0 0 0 0 0 61
62 -7.0 0 0 1 0 0 0 0 0 0 0 0 0 62
63 -12.2 0 0 0 1 0 0 0 0 0 0 0 0 63
64 -2.5 0 0 0 0 1 0 0 0 0 0 0 0 64
65 4.4 1 0 0 0 0 1 0 0 0 0 0 0 65
66 13.7 1 0 0 0 0 0 1 0 0 0 0 0 66
67 12.3 1 0 0 0 0 0 0 1 0 0 0 0 67
68 13.4 1 0 0 0 0 0 0 0 1 0 0 0 68
69 2.2 1 0 0 0 0 0 0 0 0 1 0 0 69
70 1.7 1 0 0 0 0 0 0 0 0 0 1 0 70
71 -7.2 0 0 0 0 0 0 0 0 0 0 0 1 71
72 -4.8 0 0 0 0 0 0 0 0 0 0 0 0 72
73 -2.9 0 1 0 0 0 0 0 0 0 0 0 0 73
74 -2.4 0 0 1 0 0 0 0 0 0 0 0 0 74
75 -2.5 0 0 0 1 0 0 0 0 0 0 0 0 75
76 -5.3 0 0 0 0 1 0 0 0 0 0 0 0 76
77 -7.1 0 0 0 0 0 1 0 0 0 0 0 0 77
78 -8.0 0 0 0 0 0 0 1 0 0 0 0 0 78
79 -8.9 0 0 0 0 0 0 0 1 0 0 0 0 79
80 -7.7 0 0 0 0 0 0 0 0 1 0 0 0 80
81 -1.1 0 0 0 0 0 0 0 0 0 1 0 0 81
82 4.0 1 0 0 0 0 0 0 0 0 0 1 0 82
83 9.6 1 0 0 0 0 0 0 0 0 0 0 1 83
84 10.9 1 0 0 0 0 0 0 0 0 0 0 0 84
85 13.0 1 1 0 0 0 0 0 0 0 0 0 0 85
86 14.9 1 0 1 0 0 0 0 0 0 0 0 0 86
87 20.1 1 0 0 1 0 0 0 0 0 0 0 0 87
88 10.8 1 0 0 0 1 0 0 0 0 0 0 0 88
89 11.0 1 0 0 0 0 1 0 0 0 0 0 0 89
90 3.8 1 0 0 0 0 0 1 0 0 0 0 0 90
91 10.8 1 0 0 0 0 0 0 1 0 0 0 0 91
92 7.6 1 0 0 0 0 0 0 0 1 0 0 0 92
93 10.2 1 0 0 0 0 0 0 0 0 1 0 0 93
94 2.2 1 0 0 0 0 0 0 0 0 0 1 0 94
95 -0.1 0 0 0 0 0 0 0 0 0 0 0 1 95
96 -1.7 0 0 0 0 0 0 0 0 0 0 0 0 96
97 -4.8 0 1 0 0 0 0 0 0 0 0 0 0 97
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) D M1 M2 M3 M4
-12.81751 15.26943 1.01484 1.54472 1.84650 0.24827
M5 M6 M7 M8 M9 M10
-2.08363 -2.53185 -1.61757 -2.24079 0.05717 -2.32474
M11 t
0.21072 0.09822
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-10.00979 -3.93997 -0.03981 3.26171 13.34212
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -12.81751 2.26228 -5.666 2.06e-07 ***
D 15.26943 1.17235 13.025 < 2e-16 ***
M1 1.01484 2.70323 0.375 0.708
M2 1.54472 2.78711 0.554 0.581
M3 1.84650 2.78566 0.663 0.509
M4 0.24827 2.78436 0.089 0.929
M5 -2.08363 2.78845 -0.747 0.457
M6 -2.53185 2.78726 -0.908 0.366
M7 -1.61757 2.78622 -0.581 0.563
M8 -2.24079 2.78534 -0.804 0.423
M9 0.05717 2.78016 0.021 0.984
M10 -2.32474 2.78403 -0.835 0.406
M11 0.21072 2.77955 0.076 0.940
t 0.09822 0.02062 4.764 7.98e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 5.559 on 83 degrees of freedom
Multiple R-squared: 0.729, Adjusted R-squared: 0.6865
F-statistic: 17.17 on 13 and 83 DF, p-value: < 2.2e-16
> if (n > n25) {
+ kp3 <- k + 3
+ nmkm3 <- n - k - 3
+ gqarr <- array(NA, dim=c(nmkm3-kp3+1,3))
+ numgqtests <- 0
+ numsignificant1 <- 0
+ numsignificant5 <- 0
+ numsignificant10 <- 0
+ for (mypoint in kp3:nmkm3) {
+ j <- 0
+ numgqtests <- numgqtests + 1
+ for (myalt in c('greater', 'two.sided', 'less')) {
+ j <- j + 1
+ gqarr[mypoint-kp3+1,j] <- gqtest(mylm, point=mypoint, alternative=myalt)$p.value
+ }
+ if (gqarr[mypoint-kp3+1,2] < 0.01) numsignificant1 <- numsignificant1 + 1
+ if (gqarr[mypoint-kp3+1,2] < 0.05) numsignificant5 <- numsignificant5 + 1
+ if (gqarr[mypoint-kp3+1,2] < 0.10) numsignificant10 <- numsignificant10 + 1
+ }
+ gqarr
+ }
[,1] [,2] [,3]
[1,] 0.6415744 0.7168511406 0.3584255703
[2,] 0.8281003 0.3437994690 0.1718997345
[3,] 0.9815443 0.0369113048 0.0184556524
[4,] 0.9761582 0.0476836347 0.0238418173
[5,] 0.9635626 0.0728747276 0.0364373638
[6,] 0.9682592 0.0634816030 0.0317408015
[7,] 0.9527502 0.0944996040 0.0472498020
[8,] 0.9285115 0.1429770111 0.0714885056
[9,] 0.9014289 0.1971422436 0.0985711218
[10,] 0.8624456 0.2751087018 0.1375543509
[11,] 0.8862709 0.2274582109 0.1137291055
[12,] 0.8679638 0.2640723674 0.1320361837
[13,] 0.8430814 0.3138372395 0.1569186197
[14,] 0.8087268 0.3825463716 0.1912731858
[15,] 0.8376501 0.3246997862 0.1623498931
[16,] 0.8447521 0.3104958836 0.1552479418
[17,] 0.8401666 0.3196668205 0.1598334103
[18,] 0.8131853 0.3736294530 0.1868147265
[19,] 0.9384129 0.1231741402 0.0615870701
[20,] 0.9638734 0.0722531420 0.0361265710
[21,] 0.9604517 0.0790965621 0.0395482811
[22,] 0.9658208 0.0683584035 0.0341792017
[23,] 0.9571554 0.0856891277 0.0428445638
[24,] 0.9655997 0.0688005090 0.0344002545
[25,] 0.9572766 0.0854468392 0.0427234196
[26,] 0.9482749 0.1034502266 0.0517251133
[27,] 0.9910139 0.0179721885 0.0089860943
[28,] 0.9983489 0.0033022264 0.0016511132
[29,] 0.9972485 0.0055029572 0.0027514786
[30,] 0.9984675 0.0030650350 0.0015325175
[31,] 0.9976286 0.0047428083 0.0023714042
[32,] 0.9960510 0.0078980803 0.0039490401
[33,] 0.9939656 0.0120688094 0.0060344047
[34,] 0.9925602 0.0148796297 0.0074398149
[35,] 0.9883829 0.0232341089 0.0116170544
[36,] 0.9953480 0.0093039231 0.0046519616
[37,] 0.9948050 0.0103899880 0.0051949940
[38,] 0.9914300 0.0171400082 0.0085700041
[39,] 0.9951353 0.0097294664 0.0048647332
[40,] 0.9950592 0.0098815886 0.0049407943
[41,] 0.9939542 0.0120915995 0.0060457997
[42,] 0.9933538 0.0132923236 0.0066461618
[43,] 0.9894845 0.0210309961 0.0105154981
[44,] 0.9945197 0.0109606659 0.0054803330
[45,] 0.9909327 0.0181345685 0.0090672842
[46,] 0.9869169 0.0261661591 0.0130830796
[47,] 0.9972578 0.0054844442 0.0027422221
[48,] 0.9959461 0.0081078880 0.0040539440
[49,] 0.9950575 0.0098850740 0.0049425370
[50,] 0.9980136 0.0039727421 0.0019863711
[51,] 0.9979135 0.0041729989 0.0020864995
[52,] 0.9994659 0.0010682699 0.0005341350
[53,] 0.9998646 0.0002707203 0.0001353601
[54,] 0.9996896 0.0006208669 0.0003104334
[55,] 0.9994009 0.0011982407 0.0005991203
[56,] 0.9984446 0.0031108283 0.0015554142
[57,] 0.9967264 0.0065471217 0.0032735609
[58,] 0.9922964 0.0154071179 0.0077035589
[59,] 0.9953144 0.0093711670 0.0046855835
[60,] 0.9880035 0.0239930902 0.0119965451
[61,] 0.9786682 0.0426636131 0.0213318065
[62,] 0.9574874 0.0850251955 0.0425125978
[63,] 0.9602880 0.0794240886 0.0397120443
[64,] 0.9220394 0.1559211494 0.0779605747
> postscript(file="/var/www/html/rcomp/tmp/1fkvr1227695545.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/2s8pc1227695545.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/3mlu61227695545.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/4q4s31227695545.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/5c3ea1227695545.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 = 97
Frequency = 1
1 2 3 4 5 6
8.40445689 7.57635235 7.17635235 3.77635235 -1.28996879 -3.93996879
7 8 9 10 11 12
-9.05246879 -6.42746879 -2.12364765 7.86003121 -3.04307848 -3.43057848
13 14 15 16 17 18
-3.94363560 -4.17174014 -1.27174014 -3.87174014 3.86193872 4.81193872
19 20 21 22 23 24
10.09943872 2.02443872 7.39769070 0.48136955 -5.95230930 -0.03980930
25 26 27 28 29 30
-2.45286641 -2.68097096 -9.38097096 -7.48097096 -6.54729210 -5.49729210
31 32 33 34 35 36
-10.00979210 -7.48479210 -5.98097096 4.10270790 8.46902904 7.68152904
37 38 39 40 41 42
-5.70095890 -7.02906345 1.37093655 -5.22906345 -0.29538459 2.55461541
43 44 45 46 47 48
13.34211541 12.36711541 2.27093655 6.32404624 2.09036739 0.40286739
49 50 51 52 53 54
1.58981028 4.06170574 1.16170574 9.76170574 4.49538459 0.04538459
55 56 57 58 59 60
-8.06711541 -4.84211541 3.26170574 -2.12404624 -1.95772510 -7.44522510
61 62 63 64 65 66
1.71114862 -1.81695592 -7.41695592 3.78304408 -2.35270790 7.29729210
67 68 69 70 71 72
4.88479210 6.50979210 -7.08638675 -5.30270790 -1.56695592 0.94554408
73 74 75 76 77 78
1.73248697 1.60438243 1.10438243 -0.19561757 0.23806128 -0.31193872
79 80 81 82 83 84
-2.22443872 -0.49943872 3.70438243 -4.18136955 -1.21504841 0.19745159
85 86 87 88 89 90
1.18439448 2.45628994 7.25628994 -0.54371006 1.88996879 -4.96003121
91 92 93 94 95 96
1.02746879 -1.64753121 -1.44371006 -7.16003121 3.17572077 1.68822077
97
-2.52483634
> postscript(file="/var/www/html/rcomp/tmp/6xtt51227695545.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 = 97
Frequency = 1
lag(myerror, k = 1) myerror
0 8.40445689 NA
1 7.57635235 8.40445689
2 7.17635235 7.57635235
3 3.77635235 7.17635235
4 -1.28996879 3.77635235
5 -3.93996879 -1.28996879
6 -9.05246879 -3.93996879
7 -6.42746879 -9.05246879
8 -2.12364765 -6.42746879
9 7.86003121 -2.12364765
10 -3.04307848 7.86003121
11 -3.43057848 -3.04307848
12 -3.94363560 -3.43057848
13 -4.17174014 -3.94363560
14 -1.27174014 -4.17174014
15 -3.87174014 -1.27174014
16 3.86193872 -3.87174014
17 4.81193872 3.86193872
18 10.09943872 4.81193872
19 2.02443872 10.09943872
20 7.39769070 2.02443872
21 0.48136955 7.39769070
22 -5.95230930 0.48136955
23 -0.03980930 -5.95230930
24 -2.45286641 -0.03980930
25 -2.68097096 -2.45286641
26 -9.38097096 -2.68097096
27 -7.48097096 -9.38097096
28 -6.54729210 -7.48097096
29 -5.49729210 -6.54729210
30 -10.00979210 -5.49729210
31 -7.48479210 -10.00979210
32 -5.98097096 -7.48479210
33 4.10270790 -5.98097096
34 8.46902904 4.10270790
35 7.68152904 8.46902904
36 -5.70095890 7.68152904
37 -7.02906345 -5.70095890
38 1.37093655 -7.02906345
39 -5.22906345 1.37093655
40 -0.29538459 -5.22906345
41 2.55461541 -0.29538459
42 13.34211541 2.55461541
43 12.36711541 13.34211541
44 2.27093655 12.36711541
45 6.32404624 2.27093655
46 2.09036739 6.32404624
47 0.40286739 2.09036739
48 1.58981028 0.40286739
49 4.06170574 1.58981028
50 1.16170574 4.06170574
51 9.76170574 1.16170574
52 4.49538459 9.76170574
53 0.04538459 4.49538459
54 -8.06711541 0.04538459
55 -4.84211541 -8.06711541
56 3.26170574 -4.84211541
57 -2.12404624 3.26170574
58 -1.95772510 -2.12404624
59 -7.44522510 -1.95772510
60 1.71114862 -7.44522510
61 -1.81695592 1.71114862
62 -7.41695592 -1.81695592
63 3.78304408 -7.41695592
64 -2.35270790 3.78304408
65 7.29729210 -2.35270790
66 4.88479210 7.29729210
67 6.50979210 4.88479210
68 -7.08638675 6.50979210
69 -5.30270790 -7.08638675
70 -1.56695592 -5.30270790
71 0.94554408 -1.56695592
72 1.73248697 0.94554408
73 1.60438243 1.73248697
74 1.10438243 1.60438243
75 -0.19561757 1.10438243
76 0.23806128 -0.19561757
77 -0.31193872 0.23806128
78 -2.22443872 -0.31193872
79 -0.49943872 -2.22443872
80 3.70438243 -0.49943872
81 -4.18136955 3.70438243
82 -1.21504841 -4.18136955
83 0.19745159 -1.21504841
84 1.18439448 0.19745159
85 2.45628994 1.18439448
86 7.25628994 2.45628994
87 -0.54371006 7.25628994
88 1.88996879 -0.54371006
89 -4.96003121 1.88996879
90 1.02746879 -4.96003121
91 -1.64753121 1.02746879
92 -1.44371006 -1.64753121
93 -7.16003121 -1.44371006
94 3.17572077 -7.16003121
95 1.68822077 3.17572077
96 -2.52483634 1.68822077
97 NA -2.52483634
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 7.57635235 8.40445689
[2,] 7.17635235 7.57635235
[3,] 3.77635235 7.17635235
[4,] -1.28996879 3.77635235
[5,] -3.93996879 -1.28996879
[6,] -9.05246879 -3.93996879
[7,] -6.42746879 -9.05246879
[8,] -2.12364765 -6.42746879
[9,] 7.86003121 -2.12364765
[10,] -3.04307848 7.86003121
[11,] -3.43057848 -3.04307848
[12,] -3.94363560 -3.43057848
[13,] -4.17174014 -3.94363560
[14,] -1.27174014 -4.17174014
[15,] -3.87174014 -1.27174014
[16,] 3.86193872 -3.87174014
[17,] 4.81193872 3.86193872
[18,] 10.09943872 4.81193872
[19,] 2.02443872 10.09943872
[20,] 7.39769070 2.02443872
[21,] 0.48136955 7.39769070
[22,] -5.95230930 0.48136955
[23,] -0.03980930 -5.95230930
[24,] -2.45286641 -0.03980930
[25,] -2.68097096 -2.45286641
[26,] -9.38097096 -2.68097096
[27,] -7.48097096 -9.38097096
[28,] -6.54729210 -7.48097096
[29,] -5.49729210 -6.54729210
[30,] -10.00979210 -5.49729210
[31,] -7.48479210 -10.00979210
[32,] -5.98097096 -7.48479210
[33,] 4.10270790 -5.98097096
[34,] 8.46902904 4.10270790
[35,] 7.68152904 8.46902904
[36,] -5.70095890 7.68152904
[37,] -7.02906345 -5.70095890
[38,] 1.37093655 -7.02906345
[39,] -5.22906345 1.37093655
[40,] -0.29538459 -5.22906345
[41,] 2.55461541 -0.29538459
[42,] 13.34211541 2.55461541
[43,] 12.36711541 13.34211541
[44,] 2.27093655 12.36711541
[45,] 6.32404624 2.27093655
[46,] 2.09036739 6.32404624
[47,] 0.40286739 2.09036739
[48,] 1.58981028 0.40286739
[49,] 4.06170574 1.58981028
[50,] 1.16170574 4.06170574
[51,] 9.76170574 1.16170574
[52,] 4.49538459 9.76170574
[53,] 0.04538459 4.49538459
[54,] -8.06711541 0.04538459
[55,] -4.84211541 -8.06711541
[56,] 3.26170574 -4.84211541
[57,] -2.12404624 3.26170574
[58,] -1.95772510 -2.12404624
[59,] -7.44522510 -1.95772510
[60,] 1.71114862 -7.44522510
[61,] -1.81695592 1.71114862
[62,] -7.41695592 -1.81695592
[63,] 3.78304408 -7.41695592
[64,] -2.35270790 3.78304408
[65,] 7.29729210 -2.35270790
[66,] 4.88479210 7.29729210
[67,] 6.50979210 4.88479210
[68,] -7.08638675 6.50979210
[69,] -5.30270790 -7.08638675
[70,] -1.56695592 -5.30270790
[71,] 0.94554408 -1.56695592
[72,] 1.73248697 0.94554408
[73,] 1.60438243 1.73248697
[74,] 1.10438243 1.60438243
[75,] -0.19561757 1.10438243
[76,] 0.23806128 -0.19561757
[77,] -0.31193872 0.23806128
[78,] -2.22443872 -0.31193872
[79,] -0.49943872 -2.22443872
[80,] 3.70438243 -0.49943872
[81,] -4.18136955 3.70438243
[82,] -1.21504841 -4.18136955
[83,] 0.19745159 -1.21504841
[84,] 1.18439448 0.19745159
[85,] 2.45628994 1.18439448
[86,] 7.25628994 2.45628994
[87,] -0.54371006 7.25628994
[88,] 1.88996879 -0.54371006
[89,] -4.96003121 1.88996879
[90,] 1.02746879 -4.96003121
[91,] -1.64753121 1.02746879
[92,] -1.44371006 -1.64753121
[93,] -7.16003121 -1.44371006
[94,] 3.17572077 -7.16003121
[95,] 1.68822077 3.17572077
[96,] -2.52483634 1.68822077
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 7.57635235 8.40445689
2 7.17635235 7.57635235
3 3.77635235 7.17635235
4 -1.28996879 3.77635235
5 -3.93996879 -1.28996879
6 -9.05246879 -3.93996879
7 -6.42746879 -9.05246879
8 -2.12364765 -6.42746879
9 7.86003121 -2.12364765
10 -3.04307848 7.86003121
11 -3.43057848 -3.04307848
12 -3.94363560 -3.43057848
13 -4.17174014 -3.94363560
14 -1.27174014 -4.17174014
15 -3.87174014 -1.27174014
16 3.86193872 -3.87174014
17 4.81193872 3.86193872
18 10.09943872 4.81193872
19 2.02443872 10.09943872
20 7.39769070 2.02443872
21 0.48136955 7.39769070
22 -5.95230930 0.48136955
23 -0.03980930 -5.95230930
24 -2.45286641 -0.03980930
25 -2.68097096 -2.45286641
26 -9.38097096 -2.68097096
27 -7.48097096 -9.38097096
28 -6.54729210 -7.48097096
29 -5.49729210 -6.54729210
30 -10.00979210 -5.49729210
31 -7.48479210 -10.00979210
32 -5.98097096 -7.48479210
33 4.10270790 -5.98097096
34 8.46902904 4.10270790
35 7.68152904 8.46902904
36 -5.70095890 7.68152904
37 -7.02906345 -5.70095890
38 1.37093655 -7.02906345
39 -5.22906345 1.37093655
40 -0.29538459 -5.22906345
41 2.55461541 -0.29538459
42 13.34211541 2.55461541
43 12.36711541 13.34211541
44 2.27093655 12.36711541
45 6.32404624 2.27093655
46 2.09036739 6.32404624
47 0.40286739 2.09036739
48 1.58981028 0.40286739
49 4.06170574 1.58981028
50 1.16170574 4.06170574
51 9.76170574 1.16170574
52 4.49538459 9.76170574
53 0.04538459 4.49538459
54 -8.06711541 0.04538459
55 -4.84211541 -8.06711541
56 3.26170574 -4.84211541
57 -2.12404624 3.26170574
58 -1.95772510 -2.12404624
59 -7.44522510 -1.95772510
60 1.71114862 -7.44522510
61 -1.81695592 1.71114862
62 -7.41695592 -1.81695592
63 3.78304408 -7.41695592
64 -2.35270790 3.78304408
65 7.29729210 -2.35270790
66 4.88479210 7.29729210
67 6.50979210 4.88479210
68 -7.08638675 6.50979210
69 -5.30270790 -7.08638675
70 -1.56695592 -5.30270790
71 0.94554408 -1.56695592
72 1.73248697 0.94554408
73 1.60438243 1.73248697
74 1.10438243 1.60438243
75 -0.19561757 1.10438243
76 0.23806128 -0.19561757
77 -0.31193872 0.23806128
78 -2.22443872 -0.31193872
79 -0.49943872 -2.22443872
80 3.70438243 -0.49943872
81 -4.18136955 3.70438243
82 -1.21504841 -4.18136955
83 0.19745159 -1.21504841
84 1.18439448 0.19745159
85 2.45628994 1.18439448
86 7.25628994 2.45628994
87 -0.54371006 7.25628994
88 1.88996879 -0.54371006
89 -4.96003121 1.88996879
90 1.02746879 -4.96003121
91 -1.64753121 1.02746879
92 -1.44371006 -1.64753121
93 -7.16003121 -1.44371006
94 3.17572077 -7.16003121
95 1.68822077 3.17572077
96 -2.52483634 1.68822077
> 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/7r3px1227695546.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/8o01j1227695546.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/9gulp1227695546.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/101sey1227695546.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/11pb881227695546.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/125ql01227695546.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/139g2n1227695546.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/14p8ak1227695546.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/15qngx1227695546.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/16tax61227695546.tab")
+ }
>
> system("convert tmp/1fkvr1227695545.ps tmp/1fkvr1227695545.png")
> system("convert tmp/2s8pc1227695545.ps tmp/2s8pc1227695545.png")
> system("convert tmp/3mlu61227695545.ps tmp/3mlu61227695545.png")
> system("convert tmp/4q4s31227695545.ps tmp/4q4s31227695545.png")
> system("convert tmp/5c3ea1227695545.ps tmp/5c3ea1227695545.png")
> system("convert tmp/6xtt51227695545.ps tmp/6xtt51227695545.png")
> system("convert tmp/7r3px1227695546.ps tmp/7r3px1227695546.png")
> system("convert tmp/8o01j1227695546.ps tmp/8o01j1227695546.png")
> system("convert tmp/9gulp1227695546.ps tmp/9gulp1227695546.png")
> system("convert tmp/101sey1227695546.ps tmp/101sey1227695546.png")
>
>
> proc.time()
user system elapsed
5.833 2.779 6.242