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.
Natural language support but running in an English locale
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(376.974,0,377.632,0,378.205,0,370.861,0,369.167,0,371.551,0,382.842,0,381.903,0,384.502,0,392.058,0,384.359,0,388.884,0,386.586,0,387.495,0,385.705,0,378.67,0,377.367,0,376.911,0,389.827,0,387.82,0,387.267,0,380.575,0,372.402,0,376.74,0,377.795,0,376.126,0,370.804,0,367.98,0,367.866,0,366.121,0,379.421,0,378.519,0,372.423,0,355.072,0,344.693,0,342.892,0,344.178,0,337.606,0,327.103,0,323.953,0,316.532,0,306.307,0,327.225,0,329.573,0,313.761,0,307.836,0,300.074,0,304.198,0,306.122,0,300.414,0,292.133,0,290.616,0,280.244,1,285.179,1,305.486,1,305.957,1,293.886,1,289.441,1,288.776,1,299.149,1,306.532,1,309.914,1,313.468,1,314.901,1,309.16,1,316.15,1,336.544,1,339.196,1,326.738,1,320.838,1,318.62,1,331.533,1,335.378,1),dim=c(2,73),dimnames=list(c('Maandelijksewerkloosheid','x'),1:73))
> y <- array(NA,dim=c(2,73),dimnames=list(c('Maandelijksewerkloosheid','x'),1:73))
> 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'
> #'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
Maandelijksewerkloosheid x M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 376.974 0 1 0 0 0 0 0 0 0 0 0 0
2 377.632 0 0 1 0 0 0 0 0 0 0 0 0
3 378.205 0 0 0 1 0 0 0 0 0 0 0 0
4 370.861 0 0 0 0 1 0 0 0 0 0 0 0
5 369.167 0 0 0 0 0 1 0 0 0 0 0 0
6 371.551 0 0 0 0 0 0 1 0 0 0 0 0
7 382.842 0 0 0 0 0 0 0 1 0 0 0 0
8 381.903 0 0 0 0 0 0 0 0 1 0 0 0
9 384.502 0 0 0 0 0 0 0 0 0 1 0 0
10 392.058 0 0 0 0 0 0 0 0 0 0 1 0
11 384.359 0 0 0 0 0 0 0 0 0 0 0 1
12 388.884 0 0 0 0 0 0 0 0 0 0 0 0
13 386.586 0 1 0 0 0 0 0 0 0 0 0 0
14 387.495 0 0 1 0 0 0 0 0 0 0 0 0
15 385.705 0 0 0 1 0 0 0 0 0 0 0 0
16 378.670 0 0 0 0 1 0 0 0 0 0 0 0
17 377.367 0 0 0 0 0 1 0 0 0 0 0 0
18 376.911 0 0 0 0 0 0 1 0 0 0 0 0
19 389.827 0 0 0 0 0 0 0 1 0 0 0 0
20 387.820 0 0 0 0 0 0 0 0 1 0 0 0
21 387.267 0 0 0 0 0 0 0 0 0 1 0 0
22 380.575 0 0 0 0 0 0 0 0 0 0 1 0
23 372.402 0 0 0 0 0 0 0 0 0 0 0 1
24 376.740 0 0 0 0 0 0 0 0 0 0 0 0
25 377.795 0 1 0 0 0 0 0 0 0 0 0 0
26 376.126 0 0 1 0 0 0 0 0 0 0 0 0
27 370.804 0 0 0 1 0 0 0 0 0 0 0 0
28 367.980 0 0 0 0 1 0 0 0 0 0 0 0
29 367.866 0 0 0 0 0 1 0 0 0 0 0 0
30 366.121 0 0 0 0 0 0 1 0 0 0 0 0
31 379.421 0 0 0 0 0 0 0 1 0 0 0 0
32 378.519 0 0 0 0 0 0 0 0 1 0 0 0
33 372.423 0 0 0 0 0 0 0 0 0 1 0 0
34 355.072 0 0 0 0 0 0 0 0 0 0 1 0
35 344.693 0 0 0 0 0 0 0 0 0 0 0 1
36 342.892 0 0 0 0 0 0 0 0 0 0 0 0
37 344.178 0 1 0 0 0 0 0 0 0 0 0 0
38 337.606 0 0 1 0 0 0 0 0 0 0 0 0
39 327.103 0 0 0 1 0 0 0 0 0 0 0 0
40 323.953 0 0 0 0 1 0 0 0 0 0 0 0
41 316.532 0 0 0 0 0 1 0 0 0 0 0 0
42 306.307 0 0 0 0 0 0 1 0 0 0 0 0
43 327.225 0 0 0 0 0 0 0 1 0 0 0 0
44 329.573 0 0 0 0 0 0 0 0 1 0 0 0
45 313.761 0 0 0 0 0 0 0 0 0 1 0 0
46 307.836 0 0 0 0 0 0 0 0 0 0 1 0
47 300.074 0 0 0 0 0 0 0 0 0 0 0 1
48 304.198 0 0 0 0 0 0 0 0 0 0 0 0
49 306.122 0 1 0 0 0 0 0 0 0 0 0 0
50 300.414 0 0 1 0 0 0 0 0 0 0 0 0
51 292.133 0 0 0 1 0 0 0 0 0 0 0 0
52 290.616 0 0 0 0 1 0 0 0 0 0 0 0
53 280.244 1 0 0 0 0 1 0 0 0 0 0 0
54 285.179 1 0 0 0 0 0 1 0 0 0 0 0
55 305.486 1 0 0 0 0 0 0 1 0 0 0 0
56 305.957 1 0 0 0 0 0 0 0 1 0 0 0
57 293.886 1 0 0 0 0 0 0 0 0 1 0 0
58 289.441 1 0 0 0 0 0 0 0 0 0 1 0
59 288.776 1 0 0 0 0 0 0 0 0 0 0 1
60 299.149 1 0 0 0 0 0 0 0 0 0 0 0
61 306.532 1 1 0 0 0 0 0 0 0 0 0 0
62 309.914 1 0 1 0 0 0 0 0 0 0 0 0
63 313.468 1 0 0 1 0 0 0 0 0 0 0 0
64 314.901 1 0 0 0 1 0 0 0 0 0 0 0
65 309.160 1 0 0 0 0 1 0 0 0 0 0 0
66 316.150 1 0 0 0 0 0 1 0 0 0 0 0
67 336.544 1 0 0 0 0 0 0 1 0 0 0 0
68 339.196 1 0 0 0 0 0 0 0 1 0 0 0
69 326.738 1 0 0 0 0 0 0 0 0 1 0 0
70 320.838 1 0 0 0 0 0 0 0 0 0 1 0
71 318.620 1 0 0 0 0 0 0 0 0 0 0 1
72 331.533 1 0 0 0 0 0 0 0 0 0 0 0
73 335.378 1 1 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) x M1 M2 M3 M4
356.3247 -47.2762 4.8349 -0.2475 -3.8757 -7.2819
M5 M6 M7 M8 M9 M10
-3.8433 -3.5295 12.9915 13.2620 5.8635 0.4040
M11
-5.7453
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-60.32 -18.47 11.83 20.51 35.33
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 356.3247 12.4955 28.516 < 2e-16 ***
x -47.2762 7.8347 -6.034 1.08e-07 ***
M1 4.8349 16.6566 0.290 0.773
M2 -0.2475 17.3303 -0.014 0.989
M3 -3.8757 17.3303 -0.224 0.824
M4 -7.2819 17.3303 -0.420 0.676
M5 -3.8433 17.2810 -0.222 0.825
M6 -3.5295 17.2810 -0.204 0.839
M7 12.9915 17.2810 0.752 0.455
M8 13.2620 17.2810 0.767 0.446
M9 5.8635 17.2810 0.339 0.736
M10 0.4040 17.2810 0.023 0.981
M11 -5.7453 17.2810 -0.332 0.741
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 29.93 on 60 degrees of freedom
Multiple R-squared: 0.3962, Adjusted R-squared: 0.2754
F-statistic: 3.281 on 12 and 60 DF, p-value: 0.001101
> 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,] 1.879197e-02 3.758394e-02 0.98120803
[2,] 5.665227e-03 1.133045e-02 0.99433477
[3,] 1.430124e-03 2.860247e-03 0.99856988
[4,] 4.024129e-04 8.048259e-04 0.99959759
[5,] 1.036926e-04 2.073852e-04 0.99989631
[6,] 2.379531e-05 4.759062e-05 0.99997620
[7,] 1.464766e-05 2.929532e-05 0.99998535
[8,] 9.532153e-06 1.906431e-05 0.99999047
[9,] 6.358481e-06 1.271696e-05 0.99999364
[10,] 2.083502e-06 4.167003e-06 0.99999792
[11,] 1.023788e-06 2.047576e-06 0.99999898
[12,] 1.076717e-06 2.153433e-06 0.99999892
[13,] 6.738512e-07 1.347702e-06 0.99999933
[14,] 4.453851e-07 8.907701e-07 0.99999955
[15,] 4.533333e-07 9.066665e-07 0.99999955
[16,] 3.914431e-07 7.828863e-07 0.99999961
[17,] 3.588522e-07 7.177044e-07 0.99999964
[18,] 1.818883e-06 3.637766e-06 0.99999818
[19,] 2.030437e-04 4.060874e-04 0.99979696
[20,] 4.111321e-03 8.222642e-03 0.99588868
[21,] 3.465663e-02 6.931325e-02 0.96534337
[22,] 1.063923e-01 2.127845e-01 0.89360774
[23,] 2.924293e-01 5.848586e-01 0.70757072
[24,] 5.445112e-01 9.109776e-01 0.45548882
[25,] 6.999528e-01 6.000943e-01 0.30004717
[26,] 8.381032e-01 3.237936e-01 0.16189678
[27,] 9.093793e-01 1.812413e-01 0.09062067
[28,] 9.306395e-01 1.387210e-01 0.06936050
[29,] 9.400272e-01 1.199455e-01 0.05997276
[30,] 9.550569e-01 8.988628e-02 0.04494314
[31,] 9.630064e-01 7.398715e-02 0.03699357
[32,] 9.633825e-01 7.323504e-02 0.03661752
[33,] 9.574221e-01 8.515572e-02 0.04257786
[34,] 9.478622e-01 1.042757e-01 0.05213784
[35,] 9.374795e-01 1.250410e-01 0.06252050
[36,] 9.199751e-01 1.600497e-01 0.08002487
[37,] 8.902202e-01 2.195597e-01 0.10977983
[38,] 8.539025e-01 2.921950e-01 0.14609750
[39,] 8.134503e-01 3.730994e-01 0.18654971
[40,] 7.629390e-01 4.741219e-01 0.23706096
[41,] 7.094292e-01 5.811416e-01 0.29057079
[42,] 6.384293e-01 7.231415e-01 0.36157074
> postscript(file="/var/www/html/freestat/rcomp/tmp/1s2kp1291052411.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/www/html/freestat/rcomp/tmp/2luja1291052411.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/www/html/freestat/rcomp/tmp/3luja1291052411.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/www/html/freestat/rcomp/tmp/4luja1291052411.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/www/html/freestat/rcomp/tmp/5elid1291052411.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 = 73
Frequency = 1
1 2 3 4 5 6 7
15.814380 21.554805 25.755972 21.818138 16.685610 18.755777 13.525777
8 9 10 11 12 13 14
12.316277 22.313777 35.329277 33.779610 32.559277 25.426380 31.417805
15 16 17 18 19 20 21
33.255972 29.627138 24.885610 24.115777 20.510777 18.233277 25.078777
22 23 24 25 26 27 28
23.846277 21.822610 20.415277 16.635380 20.048805 18.354972 18.937138
29 30 31 32 33 34 35
15.384610 13.325777 10.104777 8.932277 10.234777 -1.656723 -5.886390
36 37 38 39 40 41 42
-13.432723 -16.981620 -18.471195 -25.346028 -25.089862 -35.949390 -46.488223
43 44 45 46 47 48 49
-42.091223 -40.013723 -48.427223 -48.892723 -50.505390 -52.126723 -55.037620
50 51 52 53 54 55 56
-55.663195 -60.316028 -58.426862 -24.961220 -20.340054 -16.554054 -16.353554
57 58 59 60 61 62 63
-21.026054 -20.011554 -14.527220 -9.899554 -7.351450 1.112975 8.295141
64 65 66 67 68 69 70
13.134308 3.954780 10.630946 14.503946 16.885446 11.825946 11.385446
71 72 73
15.316780 22.484446 21.494550
> postscript(file="/var/www/html/freestat/rcomp/tmp/6elid1291052411.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 = 73
Frequency = 1
lag(myerror, k = 1) myerror
0 15.814380 NA
1 21.554805 15.814380
2 25.755972 21.554805
3 21.818138 25.755972
4 16.685610 21.818138
5 18.755777 16.685610
6 13.525777 18.755777
7 12.316277 13.525777
8 22.313777 12.316277
9 35.329277 22.313777
10 33.779610 35.329277
11 32.559277 33.779610
12 25.426380 32.559277
13 31.417805 25.426380
14 33.255972 31.417805
15 29.627138 33.255972
16 24.885610 29.627138
17 24.115777 24.885610
18 20.510777 24.115777
19 18.233277 20.510777
20 25.078777 18.233277
21 23.846277 25.078777
22 21.822610 23.846277
23 20.415277 21.822610
24 16.635380 20.415277
25 20.048805 16.635380
26 18.354972 20.048805
27 18.937138 18.354972
28 15.384610 18.937138
29 13.325777 15.384610
30 10.104777 13.325777
31 8.932277 10.104777
32 10.234777 8.932277
33 -1.656723 10.234777
34 -5.886390 -1.656723
35 -13.432723 -5.886390
36 -16.981620 -13.432723
37 -18.471195 -16.981620
38 -25.346028 -18.471195
39 -25.089862 -25.346028
40 -35.949390 -25.089862
41 -46.488223 -35.949390
42 -42.091223 -46.488223
43 -40.013723 -42.091223
44 -48.427223 -40.013723
45 -48.892723 -48.427223
46 -50.505390 -48.892723
47 -52.126723 -50.505390
48 -55.037620 -52.126723
49 -55.663195 -55.037620
50 -60.316028 -55.663195
51 -58.426862 -60.316028
52 -24.961220 -58.426862
53 -20.340054 -24.961220
54 -16.554054 -20.340054
55 -16.353554 -16.554054
56 -21.026054 -16.353554
57 -20.011554 -21.026054
58 -14.527220 -20.011554
59 -9.899554 -14.527220
60 -7.351450 -9.899554
61 1.112975 -7.351450
62 8.295141 1.112975
63 13.134308 8.295141
64 3.954780 13.134308
65 10.630946 3.954780
66 14.503946 10.630946
67 16.885446 14.503946
68 11.825946 16.885446
69 11.385446 11.825946
70 15.316780 11.385446
71 22.484446 15.316780
72 21.494550 22.484446
73 NA 21.494550
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 21.554805 15.814380
[2,] 25.755972 21.554805
[3,] 21.818138 25.755972
[4,] 16.685610 21.818138
[5,] 18.755777 16.685610
[6,] 13.525777 18.755777
[7,] 12.316277 13.525777
[8,] 22.313777 12.316277
[9,] 35.329277 22.313777
[10,] 33.779610 35.329277
[11,] 32.559277 33.779610
[12,] 25.426380 32.559277
[13,] 31.417805 25.426380
[14,] 33.255972 31.417805
[15,] 29.627138 33.255972
[16,] 24.885610 29.627138
[17,] 24.115777 24.885610
[18,] 20.510777 24.115777
[19,] 18.233277 20.510777
[20,] 25.078777 18.233277
[21,] 23.846277 25.078777
[22,] 21.822610 23.846277
[23,] 20.415277 21.822610
[24,] 16.635380 20.415277
[25,] 20.048805 16.635380
[26,] 18.354972 20.048805
[27,] 18.937138 18.354972
[28,] 15.384610 18.937138
[29,] 13.325777 15.384610
[30,] 10.104777 13.325777
[31,] 8.932277 10.104777
[32,] 10.234777 8.932277
[33,] -1.656723 10.234777
[34,] -5.886390 -1.656723
[35,] -13.432723 -5.886390
[36,] -16.981620 -13.432723
[37,] -18.471195 -16.981620
[38,] -25.346028 -18.471195
[39,] -25.089862 -25.346028
[40,] -35.949390 -25.089862
[41,] -46.488223 -35.949390
[42,] -42.091223 -46.488223
[43,] -40.013723 -42.091223
[44,] -48.427223 -40.013723
[45,] -48.892723 -48.427223
[46,] -50.505390 -48.892723
[47,] -52.126723 -50.505390
[48,] -55.037620 -52.126723
[49,] -55.663195 -55.037620
[50,] -60.316028 -55.663195
[51,] -58.426862 -60.316028
[52,] -24.961220 -58.426862
[53,] -20.340054 -24.961220
[54,] -16.554054 -20.340054
[55,] -16.353554 -16.554054
[56,] -21.026054 -16.353554
[57,] -20.011554 -21.026054
[58,] -14.527220 -20.011554
[59,] -9.899554 -14.527220
[60,] -7.351450 -9.899554
[61,] 1.112975 -7.351450
[62,] 8.295141 1.112975
[63,] 13.134308 8.295141
[64,] 3.954780 13.134308
[65,] 10.630946 3.954780
[66,] 14.503946 10.630946
[67,] 16.885446 14.503946
[68,] 11.825946 16.885446
[69,] 11.385446 11.825946
[70,] 15.316780 11.385446
[71,] 22.484446 15.316780
[72,] 21.494550 22.484446
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 21.554805 15.814380
2 25.755972 21.554805
3 21.818138 25.755972
4 16.685610 21.818138
5 18.755777 16.685610
6 13.525777 18.755777
7 12.316277 13.525777
8 22.313777 12.316277
9 35.329277 22.313777
10 33.779610 35.329277
11 32.559277 33.779610
12 25.426380 32.559277
13 31.417805 25.426380
14 33.255972 31.417805
15 29.627138 33.255972
16 24.885610 29.627138
17 24.115777 24.885610
18 20.510777 24.115777
19 18.233277 20.510777
20 25.078777 18.233277
21 23.846277 25.078777
22 21.822610 23.846277
23 20.415277 21.822610
24 16.635380 20.415277
25 20.048805 16.635380
26 18.354972 20.048805
27 18.937138 18.354972
28 15.384610 18.937138
29 13.325777 15.384610
30 10.104777 13.325777
31 8.932277 10.104777
32 10.234777 8.932277
33 -1.656723 10.234777
34 -5.886390 -1.656723
35 -13.432723 -5.886390
36 -16.981620 -13.432723
37 -18.471195 -16.981620
38 -25.346028 -18.471195
39 -25.089862 -25.346028
40 -35.949390 -25.089862
41 -46.488223 -35.949390
42 -42.091223 -46.488223
43 -40.013723 -42.091223
44 -48.427223 -40.013723
45 -48.892723 -48.427223
46 -50.505390 -48.892723
47 -52.126723 -50.505390
48 -55.037620 -52.126723
49 -55.663195 -55.037620
50 -60.316028 -55.663195
51 -58.426862 -60.316028
52 -24.961220 -58.426862
53 -20.340054 -24.961220
54 -16.554054 -20.340054
55 -16.353554 -16.554054
56 -21.026054 -16.353554
57 -20.011554 -21.026054
58 -14.527220 -20.011554
59 -9.899554 -14.527220
60 -7.351450 -9.899554
61 1.112975 -7.351450
62 8.295141 1.112975
63 13.134308 8.295141
64 3.954780 13.134308
65 10.630946 3.954780
66 14.503946 10.630946
67 16.885446 14.503946
68 11.825946 16.885446
69 11.385446 11.825946
70 15.316780 11.385446
71 22.484446 15.316780
72 21.494550 22.484446
> 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/7pcig1291052411.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/www/html/freestat/rcomp/tmp/8pcig1291052411.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/www/html/freestat/rcomp/tmp/9h3h11291052411.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/www/html/freestat/rcomp/tmp/10h3h11291052411.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/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/1165h51291052412.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/12gfh81291052412.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/135fd11291052412.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/14gpd41291052412.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/15jpba1291052412.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/16xzrj1291052412.tab")
+ }
>
> try(system("convert tmp/1s2kp1291052411.ps tmp/1s2kp1291052411.png",intern=TRUE))
character(0)
> try(system("convert tmp/2luja1291052411.ps tmp/2luja1291052411.png",intern=TRUE))
character(0)
> try(system("convert tmp/3luja1291052411.ps tmp/3luja1291052411.png",intern=TRUE))
character(0)
> try(system("convert tmp/4luja1291052411.ps tmp/4luja1291052411.png",intern=TRUE))
character(0)
> try(system("convert tmp/5elid1291052411.ps tmp/5elid1291052411.png",intern=TRUE))
character(0)
> try(system("convert tmp/6elid1291052411.ps tmp/6elid1291052411.png",intern=TRUE))
character(0)
> try(system("convert tmp/7pcig1291052411.ps tmp/7pcig1291052411.png",intern=TRUE))
character(0)
> try(system("convert tmp/8pcig1291052411.ps tmp/8pcig1291052411.png",intern=TRUE))
character(0)
> try(system("convert tmp/9h3h11291052411.ps tmp/9h3h11291052411.png",intern=TRUE))
character(0)
> try(system("convert tmp/10h3h11291052411.ps tmp/10h3h11291052411.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.961 2.485 4.399