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(106099,0,103235,0,98857,0,101107,0,102700,0,101477,0,99639,0,96622,0,94697,0,95093,0,112885,0,121162,0,113624,0,111632,0,106707,0,108827,0,108413,0,106249,0,104861,0,102382,0,100320,0,100228,0,117089,0,121523,0,114948,0,112831,0,107605,0,108928,1,101993,1,102850,1,99925,1,101536,1,99450,1,98305,1,110159,1,109483,1,106810,1,96279,1,91982,1,90276,1,90999,1,86622,1,83117,1,80367,1,77550,1,77443,1,92844,1,92175,1,84822,1,81632,1,78872,1,81485,1,80651,1,78192,1,76844,1,76335,1,71415,1,73899,1,86822,1,86371,1,83469,1,82662,1),dim=c(2,62),dimnames=list(c('werkl.man','Wetswijziging'),1:62))
> y <- array(NA,dim=c(2,62),dimnames=list(c('werkl.man','Wetswijziging'),1:62))
> 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
werkl.man Wetswijziging M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 106099 0 1 0 0 0 0 0 0 0 0 0 0 1
2 103235 0 0 1 0 0 0 0 0 0 0 0 0 2
3 98857 0 0 0 1 0 0 0 0 0 0 0 0 3
4 101107 0 0 0 0 1 0 0 0 0 0 0 0 4
5 102700 0 0 0 0 0 1 0 0 0 0 0 0 5
6 101477 0 0 0 0 0 0 1 0 0 0 0 0 6
7 99639 0 0 0 0 0 0 0 1 0 0 0 0 7
8 96622 0 0 0 0 0 0 0 0 1 0 0 0 8
9 94697 0 0 0 0 0 0 0 0 0 1 0 0 9
10 95093 0 0 0 0 0 0 0 0 0 0 1 0 10
11 112885 0 0 0 0 0 0 0 0 0 0 0 1 11
12 121162 0 0 0 0 0 0 0 0 0 0 0 0 12
13 113624 0 1 0 0 0 0 0 0 0 0 0 0 13
14 111632 0 0 1 0 0 0 0 0 0 0 0 0 14
15 106707 0 0 0 1 0 0 0 0 0 0 0 0 15
16 108827 0 0 0 0 1 0 0 0 0 0 0 0 16
17 108413 0 0 0 0 0 1 0 0 0 0 0 0 17
18 106249 0 0 0 0 0 0 1 0 0 0 0 0 18
19 104861 0 0 0 0 0 0 0 1 0 0 0 0 19
20 102382 0 0 0 0 0 0 0 0 1 0 0 0 20
21 100320 0 0 0 0 0 0 0 0 0 1 0 0 21
22 100228 0 0 0 0 0 0 0 0 0 0 1 0 22
23 117089 0 0 0 0 0 0 0 0 0 0 0 1 23
24 121523 0 0 0 0 0 0 0 0 0 0 0 0 24
25 114948 0 1 0 0 0 0 0 0 0 0 0 0 25
26 112831 0 0 1 0 0 0 0 0 0 0 0 0 26
27 107605 0 0 0 1 0 0 0 0 0 0 0 0 27
28 108928 1 0 0 0 1 0 0 0 0 0 0 0 28
29 101993 1 0 0 0 0 1 0 0 0 0 0 0 29
30 102850 1 0 0 0 0 0 1 0 0 0 0 0 30
31 99925 1 0 0 0 0 0 0 1 0 0 0 0 31
32 101536 1 0 0 0 0 0 0 0 1 0 0 0 32
33 99450 1 0 0 0 0 0 0 0 0 1 0 0 33
34 98305 1 0 0 0 0 0 0 0 0 0 1 0 34
35 110159 1 0 0 0 0 0 0 0 0 0 0 1 35
36 109483 1 0 0 0 0 0 0 0 0 0 0 0 36
37 106810 1 1 0 0 0 0 0 0 0 0 0 0 37
38 96279 1 0 1 0 0 0 0 0 0 0 0 0 38
39 91982 1 0 0 1 0 0 0 0 0 0 0 0 39
40 90276 1 0 0 0 1 0 0 0 0 0 0 0 40
41 90999 1 0 0 0 0 1 0 0 0 0 0 0 41
42 86622 1 0 0 0 0 0 1 0 0 0 0 0 42
43 83117 1 0 0 0 0 0 0 1 0 0 0 0 43
44 80367 1 0 0 0 0 0 0 0 1 0 0 0 44
45 77550 1 0 0 0 0 0 0 0 0 1 0 0 45
46 77443 1 0 0 0 0 0 0 0 0 0 1 0 46
47 92844 1 0 0 0 0 0 0 0 0 0 0 1 47
48 92175 1 0 0 0 0 0 0 0 0 0 0 0 48
49 84822 1 1 0 0 0 0 0 0 0 0 0 0 49
50 81632 1 0 1 0 0 0 0 0 0 0 0 0 50
51 78872 1 0 0 1 0 0 0 0 0 0 0 0 51
52 81485 1 0 0 0 1 0 0 0 0 0 0 0 52
53 80651 1 0 0 0 0 1 0 0 0 0 0 0 53
54 78192 1 0 0 0 0 0 1 0 0 0 0 0 54
55 76844 1 0 0 0 0 0 0 1 0 0 0 0 55
56 76335 1 0 0 0 0 0 0 0 1 0 0 0 56
57 71415 1 0 0 0 0 0 0 0 0 1 0 0 57
58 73899 1 0 0 0 0 0 0 0 0 0 1 0 58
59 86822 1 0 0 0 0 0 0 0 0 0 0 1 59
60 86371 1 0 0 0 0 0 0 0 0 0 0 0 60
61 83469 1 1 0 0 0 0 0 0 0 0 0 0 61
62 82662 1 0 1 0 0 0 0 0 0 0 0 0 62
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Wetswijziging M1 M2 M3
125874.13 -47.98 -7255.40 -10291.60 -14273.43
M4 M5 M6 M7 M8
-12396.54 -13022.65 -14348.56 -16002.06 -16883.57
M9 M10 M11 t
-19098.28 -18243.79 -2730.29 -547.29
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-11972 -4714 -1489 4994 11478
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 125874.13 3908.94 32.202 < 2e-16 ***
Wetswijziging -47.98 3784.82 -0.013 0.989939
M1 -7255.40 4436.65 -1.635 0.108521
M2 -10291.60 4433.21 -2.321 0.024553 *
M3 -14273.43 4649.58 -3.070 0.003519 **
M4 -12396.54 4700.97 -2.637 0.011236 *
M5 -13022.65 4683.31 -2.781 0.007727 **
M6 -14348.56 4667.94 -3.074 0.003480 **
M7 -16002.06 4654.90 -3.438 0.001222 **
M8 -16883.57 4644.20 -3.635 0.000675 ***
M9 -19098.28 4635.86 -4.120 0.000149 ***
M10 -18243.79 4629.90 -3.940 0.000263 ***
M11 -2730.29 4626.32 -0.590 0.557848
t -547.29 105.13 -5.206 3.99e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 7313 on 48 degrees of freedom
Multiple R-squared: 0.7467, Adjusted R-squared: 0.6781
F-statistic: 10.88 on 13 and 48 DF, p-value: 2.852e-10
> 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.0108040707 0.0216081415 0.989195929
[2,] 0.0098797567 0.0197595134 0.990120243
[3,] 0.0039603256 0.0079206513 0.996039674
[4,] 0.0012340211 0.0024680422 0.998765979
[5,] 0.0003851933 0.0007703866 0.999614807
[6,] 0.0001675225 0.0003350451 0.999832477
[7,] 0.0001137670 0.0002275341 0.999886233
[8,] 0.0016100897 0.0032201794 0.998389910
[9,] 0.0016221726 0.0032443452 0.998377827
[10,] 0.0009638223 0.0019276446 0.999036178
[11,] 0.0005961042 0.0011922083 0.999403896
[12,] 0.0002832709 0.0005665417 0.999716729
[13,] 0.0009824951 0.0019649902 0.999017505
[14,] 0.0004580950 0.0009161899 0.999541905
[15,] 0.0002420101 0.0004840203 0.999757990
[16,] 0.0002380931 0.0004761863 0.999761907
[17,] 0.0003944184 0.0007888367 0.999605582
[18,] 0.0004533278 0.0009066556 0.999546672
[19,] 0.0015611288 0.0031222576 0.998438871
[20,] 0.0645832447 0.1291664895 0.935416755
[21,] 0.3972775889 0.7945551779 0.602722411
[22,] 0.7369047390 0.5261905220 0.263095261
[23,] 0.9180764055 0.1638471890 0.081923595
[24,] 0.9672448593 0.0655102815 0.032755141
[25,] 0.9852909084 0.0294181832 0.014709092
[26,] 0.9901479015 0.0197041970 0.009852099
[27,] 0.9855490638 0.0289018723 0.014450936
[28,] 0.9671776696 0.0656446607 0.032822330
[29,] 0.9372773614 0.1254452772 0.062722639
> postscript(file="/var/www/html/rcomp/tmp/1fep11229537667.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/2k7nt1229537667.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/36n9t1229537667.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/4e6f71229537667.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/5ww7u1229537667.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 = 62
Frequency = 1
1 2 3 4 5 6
-11972.4405 -11252.9405 -11101.8190 -10181.4143 -7415.0143 -6764.8143
7 8 9 10 11 12
-6402.0143 -7990.2143 -7153.2143 -7064.4143 -4238.6143 1855.3857
13 14 15 16 17 18
2120.0738 3711.5738 3315.6952 4106.1000 4865.5000 4574.7000
19 20 21 22 23 24
5387.5000 4337.3000 5037.3000 4638.1000 6532.9000 8783.9000
25 26 27 28 29 30
10011.5881 11478.0881 10781.2095 10822.5905 5060.9905 7791.1905
31 32 33 34 35 36
7066.9905 10106.7905 10782.7905 9330.5905 6218.3905 3359.3905
37 38 39 40 41 42
8489.0786 1541.5786 1773.7000 -1261.8952 634.5048 -1869.2952
43 44 45 46 47 48
-3173.4952 -4494.6952 -4549.6952 -4963.8952 -4529.0952 -7381.0952
49 50 51 52 53 54
-6931.4071 -6537.9071 -4768.7857 -3485.3810 -3145.9810 -3731.7810
55 56 57 58 59 60
-2878.9810 -1959.1810 -4117.1810 -1940.3810 -3983.5810 -6617.5810
61 62
-1716.8929 1059.6071
> postscript(file="/var/www/html/rcomp/tmp/6i4bc1229537667.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 = 62
Frequency = 1
lag(myerror, k = 1) myerror
0 -11972.4405 NA
1 -11252.9405 -11972.4405
2 -11101.8190 -11252.9405
3 -10181.4143 -11101.8190
4 -7415.0143 -10181.4143
5 -6764.8143 -7415.0143
6 -6402.0143 -6764.8143
7 -7990.2143 -6402.0143
8 -7153.2143 -7990.2143
9 -7064.4143 -7153.2143
10 -4238.6143 -7064.4143
11 1855.3857 -4238.6143
12 2120.0738 1855.3857
13 3711.5738 2120.0738
14 3315.6952 3711.5738
15 4106.1000 3315.6952
16 4865.5000 4106.1000
17 4574.7000 4865.5000
18 5387.5000 4574.7000
19 4337.3000 5387.5000
20 5037.3000 4337.3000
21 4638.1000 5037.3000
22 6532.9000 4638.1000
23 8783.9000 6532.9000
24 10011.5881 8783.9000
25 11478.0881 10011.5881
26 10781.2095 11478.0881
27 10822.5905 10781.2095
28 5060.9905 10822.5905
29 7791.1905 5060.9905
30 7066.9905 7791.1905
31 10106.7905 7066.9905
32 10782.7905 10106.7905
33 9330.5905 10782.7905
34 6218.3905 9330.5905
35 3359.3905 6218.3905
36 8489.0786 3359.3905
37 1541.5786 8489.0786
38 1773.7000 1541.5786
39 -1261.8952 1773.7000
40 634.5048 -1261.8952
41 -1869.2952 634.5048
42 -3173.4952 -1869.2952
43 -4494.6952 -3173.4952
44 -4549.6952 -4494.6952
45 -4963.8952 -4549.6952
46 -4529.0952 -4963.8952
47 -7381.0952 -4529.0952
48 -6931.4071 -7381.0952
49 -6537.9071 -6931.4071
50 -4768.7857 -6537.9071
51 -3485.3810 -4768.7857
52 -3145.9810 -3485.3810
53 -3731.7810 -3145.9810
54 -2878.9810 -3731.7810
55 -1959.1810 -2878.9810
56 -4117.1810 -1959.1810
57 -1940.3810 -4117.1810
58 -3983.5810 -1940.3810
59 -6617.5810 -3983.5810
60 -1716.8929 -6617.5810
61 1059.6071 -1716.8929
62 NA 1059.6071
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -11252.9405 -11972.4405
[2,] -11101.8190 -11252.9405
[3,] -10181.4143 -11101.8190
[4,] -7415.0143 -10181.4143
[5,] -6764.8143 -7415.0143
[6,] -6402.0143 -6764.8143
[7,] -7990.2143 -6402.0143
[8,] -7153.2143 -7990.2143
[9,] -7064.4143 -7153.2143
[10,] -4238.6143 -7064.4143
[11,] 1855.3857 -4238.6143
[12,] 2120.0738 1855.3857
[13,] 3711.5738 2120.0738
[14,] 3315.6952 3711.5738
[15,] 4106.1000 3315.6952
[16,] 4865.5000 4106.1000
[17,] 4574.7000 4865.5000
[18,] 5387.5000 4574.7000
[19,] 4337.3000 5387.5000
[20,] 5037.3000 4337.3000
[21,] 4638.1000 5037.3000
[22,] 6532.9000 4638.1000
[23,] 8783.9000 6532.9000
[24,] 10011.5881 8783.9000
[25,] 11478.0881 10011.5881
[26,] 10781.2095 11478.0881
[27,] 10822.5905 10781.2095
[28,] 5060.9905 10822.5905
[29,] 7791.1905 5060.9905
[30,] 7066.9905 7791.1905
[31,] 10106.7905 7066.9905
[32,] 10782.7905 10106.7905
[33,] 9330.5905 10782.7905
[34,] 6218.3905 9330.5905
[35,] 3359.3905 6218.3905
[36,] 8489.0786 3359.3905
[37,] 1541.5786 8489.0786
[38,] 1773.7000 1541.5786
[39,] -1261.8952 1773.7000
[40,] 634.5048 -1261.8952
[41,] -1869.2952 634.5048
[42,] -3173.4952 -1869.2952
[43,] -4494.6952 -3173.4952
[44,] -4549.6952 -4494.6952
[45,] -4963.8952 -4549.6952
[46,] -4529.0952 -4963.8952
[47,] -7381.0952 -4529.0952
[48,] -6931.4071 -7381.0952
[49,] -6537.9071 -6931.4071
[50,] -4768.7857 -6537.9071
[51,] -3485.3810 -4768.7857
[52,] -3145.9810 -3485.3810
[53,] -3731.7810 -3145.9810
[54,] -2878.9810 -3731.7810
[55,] -1959.1810 -2878.9810
[56,] -4117.1810 -1959.1810
[57,] -1940.3810 -4117.1810
[58,] -3983.5810 -1940.3810
[59,] -6617.5810 -3983.5810
[60,] -1716.8929 -6617.5810
[61,] 1059.6071 -1716.8929
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -11252.9405 -11972.4405
2 -11101.8190 -11252.9405
3 -10181.4143 -11101.8190
4 -7415.0143 -10181.4143
5 -6764.8143 -7415.0143
6 -6402.0143 -6764.8143
7 -7990.2143 -6402.0143
8 -7153.2143 -7990.2143
9 -7064.4143 -7153.2143
10 -4238.6143 -7064.4143
11 1855.3857 -4238.6143
12 2120.0738 1855.3857
13 3711.5738 2120.0738
14 3315.6952 3711.5738
15 4106.1000 3315.6952
16 4865.5000 4106.1000
17 4574.7000 4865.5000
18 5387.5000 4574.7000
19 4337.3000 5387.5000
20 5037.3000 4337.3000
21 4638.1000 5037.3000
22 6532.9000 4638.1000
23 8783.9000 6532.9000
24 10011.5881 8783.9000
25 11478.0881 10011.5881
26 10781.2095 11478.0881
27 10822.5905 10781.2095
28 5060.9905 10822.5905
29 7791.1905 5060.9905
30 7066.9905 7791.1905
31 10106.7905 7066.9905
32 10782.7905 10106.7905
33 9330.5905 10782.7905
34 6218.3905 9330.5905
35 3359.3905 6218.3905
36 8489.0786 3359.3905
37 1541.5786 8489.0786
38 1773.7000 1541.5786
39 -1261.8952 1773.7000
40 634.5048 -1261.8952
41 -1869.2952 634.5048
42 -3173.4952 -1869.2952
43 -4494.6952 -3173.4952
44 -4549.6952 -4494.6952
45 -4963.8952 -4549.6952
46 -4529.0952 -4963.8952
47 -7381.0952 -4529.0952
48 -6931.4071 -7381.0952
49 -6537.9071 -6931.4071
50 -4768.7857 -6537.9071
51 -3485.3810 -4768.7857
52 -3145.9810 -3485.3810
53 -3731.7810 -3145.9810
54 -2878.9810 -3731.7810
55 -1959.1810 -2878.9810
56 -4117.1810 -1959.1810
57 -1940.3810 -4117.1810
58 -3983.5810 -1940.3810
59 -6617.5810 -3983.5810
60 -1716.8929 -6617.5810
61 1059.6071 -1716.8929
> 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/7da4t1229537667.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/81d3q1229537667.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/9e9vp1229537667.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/10qmc11229537667.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/11mqac1229537667.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/12wwmd1229537667.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/132gom1229537667.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/14l0xe1229537667.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/15h1vj1229537667.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/16sue91229537667.tab")
+ }
>
> system("convert tmp/1fep11229537667.ps tmp/1fep11229537667.png")
> system("convert tmp/2k7nt1229537667.ps tmp/2k7nt1229537667.png")
> system("convert tmp/36n9t1229537667.ps tmp/36n9t1229537667.png")
> system("convert tmp/4e6f71229537667.ps tmp/4e6f71229537667.png")
> system("convert tmp/5ww7u1229537667.ps tmp/5ww7u1229537667.png")
> system("convert tmp/6i4bc1229537667.ps tmp/6i4bc1229537667.png")
> system("convert tmp/7da4t1229537667.ps tmp/7da4t1229537667.png")
> system("convert tmp/81d3q1229537667.ps tmp/81d3q1229537667.png")
> system("convert tmp/9e9vp1229537667.ps tmp/9e9vp1229537667.png")
> system("convert tmp/10qmc11229537667.ps tmp/10qmc11229537667.png")
>
>
> proc.time()
user system elapsed
4.954 2.713 5.332