R version 2.9.0 (2009-04-17)
Copyright (C) 2009 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> x <- array(list(28029,0,29383,0,36438,0,32034,0,22679,0,24319,0,18004,0,17537,0,20366,0,22782,0,19169,0,13807,0,29743,0,25591,0,29096,0,26482,0,22405,0,27044,0,17970,0,18730,0,19684,0,19785,0,18479,0,10698,0,31956,0,29506,0,34506,0,27165,0,26736,0,23691,0,18157,0,17328,0,18205,0,20995,0,17382,0,9367,0,31124,0,26551,0,30651,0,25859,0,25100,0,25778,0,20418,0,18688,0,20424,0,24776,0,19814,1,12738,1,31566,1,30111,1,30019,1,31934,1,25826,1,26835,1,20205,1,17789,1,20520,1,22518,1,15572,1,11509,1,25447,1,24090,1,27786,1,26195,1,20516,1,22759,1,19028,1,16971,1,20036,1,22485,1),dim=c(2,70),dimnames=list(c('inschrijvingen','dummyvariabele'),1:70))
> y <- array(NA,dim=c(2,70),dimnames=list(c('inschrijvingen','dummyvariabele'),1:70))
> 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
inschrijvingen dummyvariabele M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 28029 0 1 0 0 0 0 0 0 0 0 0 0
2 29383 0 0 1 0 0 0 0 0 0 0 0 0
3 36438 0 0 0 1 0 0 0 0 0 0 0 0
4 32034 0 0 0 0 1 0 0 0 0 0 0 0
5 22679 0 0 0 0 0 1 0 0 0 0 0 0
6 24319 0 0 0 0 0 0 1 0 0 0 0 0
7 18004 0 0 0 0 0 0 0 1 0 0 0 0
8 17537 0 0 0 0 0 0 0 0 1 0 0 0
9 20366 0 0 0 0 0 0 0 0 0 1 0 0
10 22782 0 0 0 0 0 0 0 0 0 0 1 0
11 19169 0 0 0 0 0 0 0 0 0 0 0 1
12 13807 0 0 0 0 0 0 0 0 0 0 0 0
13 29743 0 1 0 0 0 0 0 0 0 0 0 0
14 25591 0 0 1 0 0 0 0 0 0 0 0 0
15 29096 0 0 0 1 0 0 0 0 0 0 0 0
16 26482 0 0 0 0 1 0 0 0 0 0 0 0
17 22405 0 0 0 0 0 1 0 0 0 0 0 0
18 27044 0 0 0 0 0 0 1 0 0 0 0 0
19 17970 0 0 0 0 0 0 0 1 0 0 0 0
20 18730 0 0 0 0 0 0 0 0 1 0 0 0
21 19684 0 0 0 0 0 0 0 0 0 1 0 0
22 19785 0 0 0 0 0 0 0 0 0 0 1 0
23 18479 0 0 0 0 0 0 0 0 0 0 0 1
24 10698 0 0 0 0 0 0 0 0 0 0 0 0
25 31956 0 1 0 0 0 0 0 0 0 0 0 0
26 29506 0 0 1 0 0 0 0 0 0 0 0 0
27 34506 0 0 0 1 0 0 0 0 0 0 0 0
28 27165 0 0 0 0 1 0 0 0 0 0 0 0
29 26736 0 0 0 0 0 1 0 0 0 0 0 0
30 23691 0 0 0 0 0 0 1 0 0 0 0 0
31 18157 0 0 0 0 0 0 0 1 0 0 0 0
32 17328 0 0 0 0 0 0 0 0 1 0 0 0
33 18205 0 0 0 0 0 0 0 0 0 1 0 0
34 20995 0 0 0 0 0 0 0 0 0 0 1 0
35 17382 0 0 0 0 0 0 0 0 0 0 0 1
36 9367 0 0 0 0 0 0 0 0 0 0 0 0
37 31124 0 1 0 0 0 0 0 0 0 0 0 0
38 26551 0 0 1 0 0 0 0 0 0 0 0 0
39 30651 0 0 0 1 0 0 0 0 0 0 0 0
40 25859 0 0 0 0 1 0 0 0 0 0 0 0
41 25100 0 0 0 0 0 1 0 0 0 0 0 0
42 25778 0 0 0 0 0 0 1 0 0 0 0 0
43 20418 0 0 0 0 0 0 0 1 0 0 0 0
44 18688 0 0 0 0 0 0 0 0 1 0 0 0
45 20424 0 0 0 0 0 0 0 0 0 1 0 0
46 24776 0 0 0 0 0 0 0 0 0 0 1 0
47 19814 1 0 0 0 0 0 0 0 0 0 0 1
48 12738 1 0 0 0 0 0 0 0 0 0 0 0
49 31566 1 1 0 0 0 0 0 0 0 0 0 0
50 30111 1 0 1 0 0 0 0 0 0 0 0 0
51 30019 1 0 0 1 0 0 0 0 0 0 0 0
52 31934 1 0 0 0 1 0 0 0 0 0 0 0
53 25826 1 0 0 0 0 1 0 0 0 0 0 0
54 26835 1 0 0 0 0 0 1 0 0 0 0 0
55 20205 1 0 0 0 0 0 0 1 0 0 0 0
56 17789 1 0 0 0 0 0 0 0 1 0 0 0
57 20520 1 0 0 0 0 0 0 0 0 1 0 0
58 22518 1 0 0 0 0 0 0 0 0 0 1 0
59 15572 1 0 0 0 0 0 0 0 0 0 0 1
60 11509 1 0 0 0 0 0 0 0 0 0 0 0
61 25447 1 1 0 0 0 0 0 0 0 0 0 0
62 24090 1 0 1 0 0 0 0 0 0 0 0 0
63 27786 1 0 0 1 0 0 0 0 0 0 0 0
64 26195 1 0 0 0 1 0 0 0 0 0 0 0
65 20516 1 0 0 0 0 1 0 0 0 0 0 0
66 22759 1 0 0 0 0 0 1 0 0 0 0 0
67 19028 1 0 0 0 0 0 0 1 0 0 0 0
68 16971 1 0 0 0 0 0 0 0 1 0 0 0
69 20036 1 0 0 0 0 0 0 0 0 1 0 0
70 22485 1 0 0 0 0 0 0 0 0 0 1 0
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) dummyvariabele M1 M2 M3
11791.5 -419.2 17992.4 15886.9 19764.3
M4 M5 M6 M7 M8
16626.4 12225.3 13419.3 7311.9 6188.8
M9 M10 M11
8220.8 10571.8 6459.4
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-3917.7 -1360.6 182.3 1359.3 4882.3
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 11791.5 962.2 12.254 < 2e-16 ***
dummyvariabele -419.2 529.2 -0.792 0.432
M1 17992.4 1271.5 14.151 < 2e-16 ***
M2 15886.9 1271.5 12.495 < 2e-16 ***
M3 19764.3 1271.5 15.544 < 2e-16 ***
M4 16626.4 1271.5 13.077 < 2e-16 ***
M5 12225.3 1271.5 9.615 1.55e-13 ***
M6 13419.3 1271.5 10.554 5.02e-15 ***
M7 7311.9 1271.5 5.751 3.67e-07 ***
M8 6188.8 1271.5 4.867 9.33e-06 ***
M9 8220.8 1271.5 6.466 2.46e-08 ***
M10 10571.8 1271.5 8.315 2.07e-11 ***
M11 6459.4 1327.5 4.866 9.38e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2099 on 57 degrees of freedom
Multiple R-squared: 0.8932, Adjusted R-squared: 0.8707
F-statistic: 39.74 on 12 and 57 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.98433552 0.03132896 0.01566448
[2,] 0.96678577 0.06642846 0.03321423
[3,] 0.95227033 0.09545935 0.04772967
[4,] 0.91677003 0.16645993 0.08322997
[5,] 0.86911243 0.26177514 0.13088757
[6,] 0.80255624 0.39488752 0.19744376
[7,] 0.79570563 0.40858875 0.20429437
[8,] 0.71825177 0.56349647 0.28174823
[9,] 0.68749727 0.62500546 0.31250273
[10,] 0.68722228 0.62555544 0.31277772
[11,] 0.65053981 0.69892038 0.34946019
[12,] 0.68717577 0.62564846 0.31282423
[13,] 0.63706731 0.72586537 0.36293269
[14,] 0.70721922 0.58556155 0.29278078
[15,] 0.66060956 0.67878088 0.33939044
[16,] 0.59325524 0.81348953 0.40674476
[17,] 0.51235128 0.97529744 0.48764872
[18,] 0.47520192 0.95040384 0.52479808
[19,] 0.42581175 0.85162349 0.57418825
[20,] 0.35802158 0.71604316 0.64197842
[21,] 0.38527771 0.77055543 0.61472229
[22,] 0.33324047 0.66648094 0.66675953
[23,] 0.27956947 0.55913893 0.72043053
[24,] 0.25129258 0.50258517 0.74870742
[25,] 0.30411472 0.60822944 0.69588528
[26,] 0.24003160 0.48006319 0.75996840
[27,] 0.17820452 0.35640904 0.82179548
[28,] 0.14283768 0.28567536 0.85716232
[29,] 0.09883061 0.19766121 0.90116939
[30,] 0.07074418 0.14148836 0.92925582
[31,] 0.06151885 0.12303770 0.93848115
[32,] 0.06010638 0.12021275 0.93989362
[33,] 0.03724084 0.07448168 0.96275916
[34,] 0.06924838 0.13849677 0.93075162
[35,] 0.14611276 0.29222552 0.85388724
[36,] 0.13965038 0.27930076 0.86034962
[37,] 0.32793855 0.65587710 0.67206145
[38,] 0.65332678 0.69334644 0.34667322
[39,] 0.95840253 0.08319495 0.04159747
> postscript(file="/var/www/html/rcomp/tmp/158g11262203222.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/29ru91262203222.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/3gvtq1262203222.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/44z5o1262203222.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/58ch21262203222.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 = 70
Frequency = 1
1 2 3 4 5 6
-1754.89831 1704.60169 4882.26836 3616.10169 -1337.73164 -891.73164
7 8 9 10 11 12
-1099.39831 -443.23164 353.76836 418.76836 918.12203 2015.52203
13 14 15 16 17 18
-40.89831 -2087.39831 -2459.73164 -1935.89831 -1611.73164 1833.26836
19 20 21 22 23 24
-1133.39831 749.76836 -328.23164 -2578.23164 228.12203 -1093.47797
25 26 27 28 29 30
2172.10169 1827.60169 2950.26836 -1252.89831 2719.26836 -1519.73164
31 32 33 34 35 36
-946.39831 -652.23164 -1807.23164 -1368.23164 -868.87797 -2424.47797
37 38 39 40 41 42
1340.10169 -1127.39831 -904.73164 -2558.89831 1083.26836 567.26836
43 44 45 46 47 48
1314.60169 707.76836 411.76836 2412.76836 1982.31695 1365.71695
49 50 51 52 53 54
2201.29661 2851.79661 -1117.53672 3935.29661 2228.46328 2043.46328
55 56 57 58 59 60
1520.79661 227.96328 926.96328 573.96328 -2259.68305 136.71695
61 62 63 64 65 66
-3917.70339 -3169.20339 -3350.53672 -1803.70339 -3081.53672 -2032.53672
67 68 69 70
343.79661 -590.03672 442.96328 540.96328
> postscript(file="/var/www/html/rcomp/tmp/6lb621262203222.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 = 70
Frequency = 1
lag(myerror, k = 1) myerror
0 -1754.89831 NA
1 1704.60169 -1754.89831
2 4882.26836 1704.60169
3 3616.10169 4882.26836
4 -1337.73164 3616.10169
5 -891.73164 -1337.73164
6 -1099.39831 -891.73164
7 -443.23164 -1099.39831
8 353.76836 -443.23164
9 418.76836 353.76836
10 918.12203 418.76836
11 2015.52203 918.12203
12 -40.89831 2015.52203
13 -2087.39831 -40.89831
14 -2459.73164 -2087.39831
15 -1935.89831 -2459.73164
16 -1611.73164 -1935.89831
17 1833.26836 -1611.73164
18 -1133.39831 1833.26836
19 749.76836 -1133.39831
20 -328.23164 749.76836
21 -2578.23164 -328.23164
22 228.12203 -2578.23164
23 -1093.47797 228.12203
24 2172.10169 -1093.47797
25 1827.60169 2172.10169
26 2950.26836 1827.60169
27 -1252.89831 2950.26836
28 2719.26836 -1252.89831
29 -1519.73164 2719.26836
30 -946.39831 -1519.73164
31 -652.23164 -946.39831
32 -1807.23164 -652.23164
33 -1368.23164 -1807.23164
34 -868.87797 -1368.23164
35 -2424.47797 -868.87797
36 1340.10169 -2424.47797
37 -1127.39831 1340.10169
38 -904.73164 -1127.39831
39 -2558.89831 -904.73164
40 1083.26836 -2558.89831
41 567.26836 1083.26836
42 1314.60169 567.26836
43 707.76836 1314.60169
44 411.76836 707.76836
45 2412.76836 411.76836
46 1982.31695 2412.76836
47 1365.71695 1982.31695
48 2201.29661 1365.71695
49 2851.79661 2201.29661
50 -1117.53672 2851.79661
51 3935.29661 -1117.53672
52 2228.46328 3935.29661
53 2043.46328 2228.46328
54 1520.79661 2043.46328
55 227.96328 1520.79661
56 926.96328 227.96328
57 573.96328 926.96328
58 -2259.68305 573.96328
59 136.71695 -2259.68305
60 -3917.70339 136.71695
61 -3169.20339 -3917.70339
62 -3350.53672 -3169.20339
63 -1803.70339 -3350.53672
64 -3081.53672 -1803.70339
65 -2032.53672 -3081.53672
66 343.79661 -2032.53672
67 -590.03672 343.79661
68 442.96328 -590.03672
69 540.96328 442.96328
70 NA 540.96328
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 1704.60169 -1754.89831
[2,] 4882.26836 1704.60169
[3,] 3616.10169 4882.26836
[4,] -1337.73164 3616.10169
[5,] -891.73164 -1337.73164
[6,] -1099.39831 -891.73164
[7,] -443.23164 -1099.39831
[8,] 353.76836 -443.23164
[9,] 418.76836 353.76836
[10,] 918.12203 418.76836
[11,] 2015.52203 918.12203
[12,] -40.89831 2015.52203
[13,] -2087.39831 -40.89831
[14,] -2459.73164 -2087.39831
[15,] -1935.89831 -2459.73164
[16,] -1611.73164 -1935.89831
[17,] 1833.26836 -1611.73164
[18,] -1133.39831 1833.26836
[19,] 749.76836 -1133.39831
[20,] -328.23164 749.76836
[21,] -2578.23164 -328.23164
[22,] 228.12203 -2578.23164
[23,] -1093.47797 228.12203
[24,] 2172.10169 -1093.47797
[25,] 1827.60169 2172.10169
[26,] 2950.26836 1827.60169
[27,] -1252.89831 2950.26836
[28,] 2719.26836 -1252.89831
[29,] -1519.73164 2719.26836
[30,] -946.39831 -1519.73164
[31,] -652.23164 -946.39831
[32,] -1807.23164 -652.23164
[33,] -1368.23164 -1807.23164
[34,] -868.87797 -1368.23164
[35,] -2424.47797 -868.87797
[36,] 1340.10169 -2424.47797
[37,] -1127.39831 1340.10169
[38,] -904.73164 -1127.39831
[39,] -2558.89831 -904.73164
[40,] 1083.26836 -2558.89831
[41,] 567.26836 1083.26836
[42,] 1314.60169 567.26836
[43,] 707.76836 1314.60169
[44,] 411.76836 707.76836
[45,] 2412.76836 411.76836
[46,] 1982.31695 2412.76836
[47,] 1365.71695 1982.31695
[48,] 2201.29661 1365.71695
[49,] 2851.79661 2201.29661
[50,] -1117.53672 2851.79661
[51,] 3935.29661 -1117.53672
[52,] 2228.46328 3935.29661
[53,] 2043.46328 2228.46328
[54,] 1520.79661 2043.46328
[55,] 227.96328 1520.79661
[56,] 926.96328 227.96328
[57,] 573.96328 926.96328
[58,] -2259.68305 573.96328
[59,] 136.71695 -2259.68305
[60,] -3917.70339 136.71695
[61,] -3169.20339 -3917.70339
[62,] -3350.53672 -3169.20339
[63,] -1803.70339 -3350.53672
[64,] -3081.53672 -1803.70339
[65,] -2032.53672 -3081.53672
[66,] 343.79661 -2032.53672
[67,] -590.03672 343.79661
[68,] 442.96328 -590.03672
[69,] 540.96328 442.96328
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 1704.60169 -1754.89831
2 4882.26836 1704.60169
3 3616.10169 4882.26836
4 -1337.73164 3616.10169
5 -891.73164 -1337.73164
6 -1099.39831 -891.73164
7 -443.23164 -1099.39831
8 353.76836 -443.23164
9 418.76836 353.76836
10 918.12203 418.76836
11 2015.52203 918.12203
12 -40.89831 2015.52203
13 -2087.39831 -40.89831
14 -2459.73164 -2087.39831
15 -1935.89831 -2459.73164
16 -1611.73164 -1935.89831
17 1833.26836 -1611.73164
18 -1133.39831 1833.26836
19 749.76836 -1133.39831
20 -328.23164 749.76836
21 -2578.23164 -328.23164
22 228.12203 -2578.23164
23 -1093.47797 228.12203
24 2172.10169 -1093.47797
25 1827.60169 2172.10169
26 2950.26836 1827.60169
27 -1252.89831 2950.26836
28 2719.26836 -1252.89831
29 -1519.73164 2719.26836
30 -946.39831 -1519.73164
31 -652.23164 -946.39831
32 -1807.23164 -652.23164
33 -1368.23164 -1807.23164
34 -868.87797 -1368.23164
35 -2424.47797 -868.87797
36 1340.10169 -2424.47797
37 -1127.39831 1340.10169
38 -904.73164 -1127.39831
39 -2558.89831 -904.73164
40 1083.26836 -2558.89831
41 567.26836 1083.26836
42 1314.60169 567.26836
43 707.76836 1314.60169
44 411.76836 707.76836
45 2412.76836 411.76836
46 1982.31695 2412.76836
47 1365.71695 1982.31695
48 2201.29661 1365.71695
49 2851.79661 2201.29661
50 -1117.53672 2851.79661
51 3935.29661 -1117.53672
52 2228.46328 3935.29661
53 2043.46328 2228.46328
54 1520.79661 2043.46328
55 227.96328 1520.79661
56 926.96328 227.96328
57 573.96328 926.96328
58 -2259.68305 573.96328
59 136.71695 -2259.68305
60 -3917.70339 136.71695
61 -3169.20339 -3917.70339
62 -3350.53672 -3169.20339
63 -1803.70339 -3350.53672
64 -3081.53672 -1803.70339
65 -2032.53672 -3081.53672
66 343.79661 -2032.53672
67 -590.03672 343.79661
68 442.96328 -590.03672
69 540.96328 442.96328
> 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/7ce3w1262203222.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/84byo1262203222.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/9vukh1262203222.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/10y0op1262203222.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/11o1lq1262203222.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/12gyp01262203222.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/1336nj1262203222.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/14wae51262203223.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/15ombq1262203223.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/16rva81262203223.tab")
+ }
>
> try(system("convert tmp/158g11262203222.ps tmp/158g11262203222.png",intern=TRUE))
character(0)
> try(system("convert tmp/29ru91262203222.ps tmp/29ru91262203222.png",intern=TRUE))
character(0)
> try(system("convert tmp/3gvtq1262203222.ps tmp/3gvtq1262203222.png",intern=TRUE))
character(0)
> try(system("convert tmp/44z5o1262203222.ps tmp/44z5o1262203222.png",intern=TRUE))
character(0)
> try(system("convert tmp/58ch21262203222.ps tmp/58ch21262203222.png",intern=TRUE))
character(0)
> try(system("convert tmp/6lb621262203222.ps tmp/6lb621262203222.png",intern=TRUE))
character(0)
> try(system("convert tmp/7ce3w1262203222.ps tmp/7ce3w1262203222.png",intern=TRUE))
character(0)
> try(system("convert tmp/84byo1262203222.ps tmp/84byo1262203222.png",intern=TRUE))
character(0)
> try(system("convert tmp/9vukh1262203222.ps tmp/9vukh1262203222.png",intern=TRUE))
character(0)
> try(system("convert tmp/10y0op1262203222.ps tmp/10y0op1262203222.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
2.554 1.589 3.242