R version 2.12.0 (2010-10-15)
Copyright (C) 2010 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: i486-pc-linux-gnu (32-bit)
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(63047
+ ,13
+ ,6
+ ,10823
+ ,66751
+ ,26
+ ,7
+ ,44480
+ ,7176
+ ,0
+ ,0
+ ,1929
+ ,78306
+ ,37
+ ,12
+ ,30032
+ ,144655
+ ,47
+ ,15
+ ,27669
+ ,269638
+ ,84
+ ,16
+ ,114967
+ ,69266
+ ,21
+ ,12
+ ,29951
+ ,83529
+ ,36
+ ,15
+ ,38824
+ ,73226
+ ,35
+ ,15
+ ,26517
+ ,178519
+ ,40
+ ,13
+ ,63570
+ ,67250
+ ,35
+ ,6
+ ,27131
+ ,102982
+ ,46
+ ,16
+ ,41061
+ ,50001
+ ,20
+ ,7
+ ,18810
+ ,91093
+ ,24
+ ,12
+ ,27582
+ ,80112
+ ,19
+ ,9
+ ,37026
+ ,72961
+ ,15
+ ,10
+ ,24252
+ ,77159
+ ,52
+ ,16
+ ,32579
+ ,15629
+ ,0
+ ,5
+ ,0
+ ,71693
+ ,38
+ ,20
+ ,29666
+ ,19920
+ ,12
+ ,7
+ ,7533
+ ,39403
+ ,10
+ ,13
+ ,11892
+ ,104383
+ ,53
+ ,13
+ ,51557
+ ,56088
+ ,4
+ ,11
+ ,5737
+ ,62006
+ ,24
+ ,9
+ ,11203
+ ,81665
+ ,39
+ ,10
+ ,28714
+ ,69451
+ ,19
+ ,7
+ ,24268
+ ,88794
+ ,23
+ ,13
+ ,30749
+ ,90642
+ ,39
+ ,15
+ ,46643
+ ,207069
+ ,38
+ ,13
+ ,64530
+ ,99340
+ ,20
+ ,7
+ ,35346
+ ,56695
+ ,20
+ ,14
+ ,5766
+ ,108143
+ ,41
+ ,11
+ ,29217
+ ,64204
+ ,29
+ ,3
+ ,15912
+ ,29101
+ ,0
+ ,8
+ ,3728
+ ,113060
+ ,31
+ ,12
+ ,37494
+ ,0
+ ,0
+ ,0
+ ,0
+ ,65773
+ ,8
+ ,12
+ ,13214
+ ,67047
+ ,35
+ ,8
+ ,19576
+ ,41953
+ ,3
+ ,20
+ ,13632
+ ,113787
+ ,47
+ ,18
+ ,67378
+ ,86584
+ ,42
+ ,9
+ ,29387
+ ,59588
+ ,11
+ ,14
+ ,15936
+ ,40064
+ ,10
+ ,7
+ ,18156
+ ,74471
+ ,26
+ ,13
+ ,23750
+ ,60437
+ ,27
+ ,11
+ ,15559
+ ,55118
+ ,1
+ ,11
+ ,21713
+ ,40295
+ ,15
+ ,14
+ ,12023
+ ,103397
+ ,32
+ ,9
+ ,23588
+ ,78982
+ ,13
+ ,12
+ ,28661
+ ,67317
+ ,25
+ ,12
+ ,16874
+ ,39887
+ ,10
+ ,17
+ ,11804
+ ,59682
+ ,24
+ ,10
+ ,12949
+ ,132740
+ ,26
+ ,11
+ ,38340
+ ,104816
+ ,29
+ ,12
+ ,36573
+ ,101395
+ ,40
+ ,17
+ ,40068
+ ,72824
+ ,22
+ ,6
+ ,25561
+ ,76018
+ ,27
+ ,8
+ ,31287
+ ,33891
+ ,8
+ ,12
+ ,8383
+ ,63694
+ ,27
+ ,13
+ ,29178
+ ,33239
+ ,0
+ ,14
+ ,1237
+ ,35093
+ ,0
+ ,17
+ ,10241
+ ,35252
+ ,17
+ ,8
+ ,8219
+ ,36977
+ ,7
+ ,9
+ ,9348
+ ,42406
+ ,18
+ ,9
+ ,25242
+ ,56353
+ ,7
+ ,9
+ ,24267
+ ,58817
+ ,24
+ ,15
+ ,25902
+ ,81051
+ ,19
+ ,16
+ ,51849
+ ,70872
+ ,39
+ ,13
+ ,29065
+ ,42372
+ ,17
+ ,12
+ ,22417)
+ ,dim=c(4
+ ,69)
+ ,dimnames=list(c('TijdInRFC'
+ ,'Blogs'
+ ,'Peerreviews'
+ ,'Compendiumtijd')
+ ,1:69))
> y <- array(NA,dim=c(4,69),dimnames=list(c('TijdInRFC','Blogs','Peerreviews','Compendiumtijd'),1:69))
> 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 = 'Do not include Seasonal 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
> 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
TijdInRFC Blogs Peerreviews Compendiumtijd
1 63047 13 6 10823
2 66751 26 7 44480
3 7176 0 0 1929
4 78306 37 12 30032
5 144655 47 15 27669
6 269638 84 16 114967
7 69266 21 12 29951
8 83529 36 15 38824
9 73226 35 15 26517
10 178519 40 13 63570
11 67250 35 6 27131
12 102982 46 16 41061
13 50001 20 7 18810
14 91093 24 12 27582
15 80112 19 9 37026
16 72961 15 10 24252
17 77159 52 16 32579
18 15629 0 5 0
19 71693 38 20 29666
20 19920 12 7 7533
21 39403 10 13 11892
22 104383 53 13 51557
23 56088 4 11 5737
24 62006 24 9 11203
25 81665 39 10 28714
26 69451 19 7 24268
27 88794 23 13 30749
28 90642 39 15 46643
29 207069 38 13 64530
30 99340 20 7 35346
31 56695 20 14 5766
32 108143 41 11 29217
33 64204 29 3 15912
34 29101 0 8 3728
35 113060 31 12 37494
36 0 0 0 0
37 65773 8 12 13214
38 67047 35 8 19576
39 41953 3 20 13632
40 113787 47 18 67378
41 86584 42 9 29387
42 59588 11 14 15936
43 40064 10 7 18156
44 74471 26 13 23750
45 60437 27 11 15559
46 55118 1 11 21713
47 40295 15 14 12023
48 103397 32 9 23588
49 78982 13 12 28661
50 67317 25 12 16874
51 39887 10 17 11804
52 59682 24 10 12949
53 132740 26 11 38340
54 104816 29 12 36573
55 101395 40 17 40068
56 72824 22 6 25561
57 76018 27 8 31287
58 33891 8 12 8383
59 63694 27 13 29178
60 33239 0 14 1237
61 35093 0 17 10241
62 35252 17 8 8219
63 36977 7 9 9348
64 42406 18 9 25242
65 56353 7 9 24267
66 58817 24 15 25902
67 81051 19 16 51849
68 70872 39 13 29065
69 42372 17 12 22417
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Blogs Peerreviews Compendiumtijd
16866.061 592.162 59.988 1.613
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-40638 -12420 -1950 9671 62866
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.687e+04 7.048e+03 2.393 0.0196 *
Blogs 5.922e+02 2.415e+02 2.452 0.0169 *
Peerreviews 5.999e+01 6.273e+02 0.096 0.9241
Compendiumtijd 1.613e+00 2.102e-01 7.671 1.12e-10 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 19610 on 65 degrees of freedom
Multiple R-squared: 0.7966, Adjusted R-squared: 0.7872
F-statistic: 84.86 on 3 and 65 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.8827200 0.234560051 0.117280025
[2,] 0.8485951 0.302809708 0.151404854
[3,] 0.8078127 0.384374667 0.192187333
[4,] 0.9789617 0.042076592 0.021038296
[5,] 0.9777209 0.044558271 0.022279135
[6,] 0.9714796 0.057040826 0.028520413
[7,] 0.9533981 0.093203805 0.046601903
[8,] 0.9433719 0.113256140 0.056628070
[9,] 0.9141575 0.171684941 0.085842471
[10,] 0.8862977 0.227404547 0.113702274
[11,] 0.9208366 0.158326894 0.079163447
[12,] 0.8888563 0.222287498 0.111143749
[13,] 0.8775203 0.244959352 0.122479676
[14,] 0.8545688 0.290862312 0.145431156
[15,] 0.8051905 0.389618920 0.194809460
[16,] 0.8401585 0.319683086 0.159841543
[17,] 0.8718515 0.256297030 0.128148515
[18,] 0.8519811 0.296037893 0.148018946
[19,] 0.8060117 0.387976648 0.193988324
[20,] 0.7517290 0.496541962 0.248270981
[21,] 0.6954885 0.609023034 0.304511517
[22,] 0.7367805 0.526439004 0.263219502
[23,] 0.9901229 0.019754144 0.009877072
[24,] 0.9892945 0.021410910 0.010705455
[25,] 0.9867104 0.026579104 0.013289552
[26,] 0.9881075 0.023784995 0.011892498
[27,] 0.9823632 0.035273512 0.017636756
[28,] 0.9728196 0.054360718 0.027180359
[29,] 0.9754796 0.049040752 0.024520376
[30,] 0.9805189 0.038962130 0.019481065
[31,] 0.9815382 0.036923530 0.018461765
[32,] 0.9723135 0.055372957 0.027686478
[33,] 0.9616194 0.076761198 0.038380599
[34,] 0.9815038 0.036992456 0.018496228
[35,] 0.9714270 0.057146021 0.028573011
[36,] 0.9620012 0.075997520 0.037998760
[37,] 0.9541355 0.091728926 0.045864463
[38,] 0.9328291 0.134341828 0.067170914
[39,] 0.9025603 0.194879383 0.097439691
[40,] 0.8668556 0.266288757 0.133144378
[41,] 0.8221648 0.355670318 0.177835159
[42,] 0.8811811 0.237637863 0.118818932
[43,] 0.8541302 0.291739546 0.145869773
[44,] 0.8111828 0.377634318 0.188817159
[45,] 0.7463651 0.507269829 0.253634915
[46,] 0.6817081 0.636583851 0.318291926
[47,] 0.9775208 0.044958329 0.022479164
[48,] 0.9950278 0.009944326 0.004972163
[49,] 0.9977428 0.004514478 0.002257239
[50,] 0.9981347 0.003730587 0.001865294
[51,] 0.9989469 0.002106221 0.001053110
[52,] 0.9967426 0.006514779 0.003257390
[53,] 0.9909162 0.018167610 0.009083805
[54,] 0.9794768 0.041046477 0.020523239
[55,] 0.9433720 0.113255961 0.056627981
[56,] 0.8571219 0.285756220 0.142878110
> postscript(file="/var/www/rcomp/tmp/1spch1322155631.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/rcomp/tmp/2meoy1322155631.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/rcomp/tmp/3c6ck1322155631.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/rcomp/tmp/4yysq1322155631.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/rcomp/tmp/5otoq1322155631.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 = 69
Frequency = 1
1 2 3 4 5 6
20670.7938 -37655.2777 -12800.5774 -9616.5759 54441.1794 16686.0030
7 8 9 10 11 12
-9051.3676 -18158.4941 -8024.2717 34679.8625 -14450.4553 -8294.2711
13 14 15 16 17 18
-9459.3820 14819.1626 -8249.5295 7506.2299 -23993.0039 -1537.0016
19 20 21 22 23 24
-16711.4670 -16618.9009 -3340.4003 -27783.2619 26942.4980 12323.2944
25 26 27 28 29 30
-5196.6469 1781.7451 7945.5434 -25430.1327 62866.1852 13215.2875
31 32 33 34 35 36
17848.1756 19225.9522 4327.1404 5743.6273 16657.9087 -16866.0614
37 38 39 40 41 42
22142.1848 -2590.9800 129.0657 -40637.6204 -3079.3592 9671.4919
43 44 45 46 47 48
-12420.1835 3131.9571 1833.7736 1987.6553 -5680.4373 29006.1613
49 50 51 52 53 54
7482.0577 7717.6701 -2954.4524 7123.8781 37994.5314 11083.3474
55 56 57 58 59 60
-4787.0714 1353.2816 -7766.6912 -1949.8188 -16989.8651 13538.4411
61 62 63 64 65 66
693.5103 -5413.8764 352.2433 -26361.6455 -4328.6717 -14927.7985
67 68 69
-31632.5607 -16735.5992 -21428.1298
> postscript(file="/var/www/rcomp/tmp/66m0w1322155631.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 = 69
Frequency = 1
lag(myerror, k = 1) myerror
0 20670.7938 NA
1 -37655.2777 20670.7938
2 -12800.5774 -37655.2777
3 -9616.5759 -12800.5774
4 54441.1794 -9616.5759
5 16686.0030 54441.1794
6 -9051.3676 16686.0030
7 -18158.4941 -9051.3676
8 -8024.2717 -18158.4941
9 34679.8625 -8024.2717
10 -14450.4553 34679.8625
11 -8294.2711 -14450.4553
12 -9459.3820 -8294.2711
13 14819.1626 -9459.3820
14 -8249.5295 14819.1626
15 7506.2299 -8249.5295
16 -23993.0039 7506.2299
17 -1537.0016 -23993.0039
18 -16711.4670 -1537.0016
19 -16618.9009 -16711.4670
20 -3340.4003 -16618.9009
21 -27783.2619 -3340.4003
22 26942.4980 -27783.2619
23 12323.2944 26942.4980
24 -5196.6469 12323.2944
25 1781.7451 -5196.6469
26 7945.5434 1781.7451
27 -25430.1327 7945.5434
28 62866.1852 -25430.1327
29 13215.2875 62866.1852
30 17848.1756 13215.2875
31 19225.9522 17848.1756
32 4327.1404 19225.9522
33 5743.6273 4327.1404
34 16657.9087 5743.6273
35 -16866.0614 16657.9087
36 22142.1848 -16866.0614
37 -2590.9800 22142.1848
38 129.0657 -2590.9800
39 -40637.6204 129.0657
40 -3079.3592 -40637.6204
41 9671.4919 -3079.3592
42 -12420.1835 9671.4919
43 3131.9571 -12420.1835
44 1833.7736 3131.9571
45 1987.6553 1833.7736
46 -5680.4373 1987.6553
47 29006.1613 -5680.4373
48 7482.0577 29006.1613
49 7717.6701 7482.0577
50 -2954.4524 7717.6701
51 7123.8781 -2954.4524
52 37994.5314 7123.8781
53 11083.3474 37994.5314
54 -4787.0714 11083.3474
55 1353.2816 -4787.0714
56 -7766.6912 1353.2816
57 -1949.8188 -7766.6912
58 -16989.8651 -1949.8188
59 13538.4411 -16989.8651
60 693.5103 13538.4411
61 -5413.8764 693.5103
62 352.2433 -5413.8764
63 -26361.6455 352.2433
64 -4328.6717 -26361.6455
65 -14927.7985 -4328.6717
66 -31632.5607 -14927.7985
67 -16735.5992 -31632.5607
68 -21428.1298 -16735.5992
69 NA -21428.1298
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -37655.2777 20670.7938
[2,] -12800.5774 -37655.2777
[3,] -9616.5759 -12800.5774
[4,] 54441.1794 -9616.5759
[5,] 16686.0030 54441.1794
[6,] -9051.3676 16686.0030
[7,] -18158.4941 -9051.3676
[8,] -8024.2717 -18158.4941
[9,] 34679.8625 -8024.2717
[10,] -14450.4553 34679.8625
[11,] -8294.2711 -14450.4553
[12,] -9459.3820 -8294.2711
[13,] 14819.1626 -9459.3820
[14,] -8249.5295 14819.1626
[15,] 7506.2299 -8249.5295
[16,] -23993.0039 7506.2299
[17,] -1537.0016 -23993.0039
[18,] -16711.4670 -1537.0016
[19,] -16618.9009 -16711.4670
[20,] -3340.4003 -16618.9009
[21,] -27783.2619 -3340.4003
[22,] 26942.4980 -27783.2619
[23,] 12323.2944 26942.4980
[24,] -5196.6469 12323.2944
[25,] 1781.7451 -5196.6469
[26,] 7945.5434 1781.7451
[27,] -25430.1327 7945.5434
[28,] 62866.1852 -25430.1327
[29,] 13215.2875 62866.1852
[30,] 17848.1756 13215.2875
[31,] 19225.9522 17848.1756
[32,] 4327.1404 19225.9522
[33,] 5743.6273 4327.1404
[34,] 16657.9087 5743.6273
[35,] -16866.0614 16657.9087
[36,] 22142.1848 -16866.0614
[37,] -2590.9800 22142.1848
[38,] 129.0657 -2590.9800
[39,] -40637.6204 129.0657
[40,] -3079.3592 -40637.6204
[41,] 9671.4919 -3079.3592
[42,] -12420.1835 9671.4919
[43,] 3131.9571 -12420.1835
[44,] 1833.7736 3131.9571
[45,] 1987.6553 1833.7736
[46,] -5680.4373 1987.6553
[47,] 29006.1613 -5680.4373
[48,] 7482.0577 29006.1613
[49,] 7717.6701 7482.0577
[50,] -2954.4524 7717.6701
[51,] 7123.8781 -2954.4524
[52,] 37994.5314 7123.8781
[53,] 11083.3474 37994.5314
[54,] -4787.0714 11083.3474
[55,] 1353.2816 -4787.0714
[56,] -7766.6912 1353.2816
[57,] -1949.8188 -7766.6912
[58,] -16989.8651 -1949.8188
[59,] 13538.4411 -16989.8651
[60,] 693.5103 13538.4411
[61,] -5413.8764 693.5103
[62,] 352.2433 -5413.8764
[63,] -26361.6455 352.2433
[64,] -4328.6717 -26361.6455
[65,] -14927.7985 -4328.6717
[66,] -31632.5607 -14927.7985
[67,] -16735.5992 -31632.5607
[68,] -21428.1298 -16735.5992
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -37655.2777 20670.7938
2 -12800.5774 -37655.2777
3 -9616.5759 -12800.5774
4 54441.1794 -9616.5759
5 16686.0030 54441.1794
6 -9051.3676 16686.0030
7 -18158.4941 -9051.3676
8 -8024.2717 -18158.4941
9 34679.8625 -8024.2717
10 -14450.4553 34679.8625
11 -8294.2711 -14450.4553
12 -9459.3820 -8294.2711
13 14819.1626 -9459.3820
14 -8249.5295 14819.1626
15 7506.2299 -8249.5295
16 -23993.0039 7506.2299
17 -1537.0016 -23993.0039
18 -16711.4670 -1537.0016
19 -16618.9009 -16711.4670
20 -3340.4003 -16618.9009
21 -27783.2619 -3340.4003
22 26942.4980 -27783.2619
23 12323.2944 26942.4980
24 -5196.6469 12323.2944
25 1781.7451 -5196.6469
26 7945.5434 1781.7451
27 -25430.1327 7945.5434
28 62866.1852 -25430.1327
29 13215.2875 62866.1852
30 17848.1756 13215.2875
31 19225.9522 17848.1756
32 4327.1404 19225.9522
33 5743.6273 4327.1404
34 16657.9087 5743.6273
35 -16866.0614 16657.9087
36 22142.1848 -16866.0614
37 -2590.9800 22142.1848
38 129.0657 -2590.9800
39 -40637.6204 129.0657
40 -3079.3592 -40637.6204
41 9671.4919 -3079.3592
42 -12420.1835 9671.4919
43 3131.9571 -12420.1835
44 1833.7736 3131.9571
45 1987.6553 1833.7736
46 -5680.4373 1987.6553
47 29006.1613 -5680.4373
48 7482.0577 29006.1613
49 7717.6701 7482.0577
50 -2954.4524 7717.6701
51 7123.8781 -2954.4524
52 37994.5314 7123.8781
53 11083.3474 37994.5314
54 -4787.0714 11083.3474
55 1353.2816 -4787.0714
56 -7766.6912 1353.2816
57 -1949.8188 -7766.6912
58 -16989.8651 -1949.8188
59 13538.4411 -16989.8651
60 693.5103 13538.4411
61 -5413.8764 693.5103
62 352.2433 -5413.8764
63 -26361.6455 352.2433
64 -4328.6717 -26361.6455
65 -14927.7985 -4328.6717
66 -31632.5607 -14927.7985
67 -16735.5992 -31632.5607
68 -21428.1298 -16735.5992
> 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/rcomp/tmp/77j6x1322155631.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/rcomp/tmp/8p3vh1322155631.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/rcomp/tmp/9s7zf1322155631.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/rcomp/tmp/10moey1322155631.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/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/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/rcomp/tmp/11vqit1322155631.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/rcomp/tmp/12pbjd1322155631.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/rcomp/tmp/13584b1322155631.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/rcomp/tmp/14u4811322155631.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/rcomp/tmp/15ofye1322155631.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/rcomp/tmp/16iuwb1322155631.tab")
+ }
>
> try(system("convert tmp/1spch1322155631.ps tmp/1spch1322155631.png",intern=TRUE))
character(0)
> try(system("convert tmp/2meoy1322155631.ps tmp/2meoy1322155631.png",intern=TRUE))
character(0)
> try(system("convert tmp/3c6ck1322155631.ps tmp/3c6ck1322155631.png",intern=TRUE))
character(0)
> try(system("convert tmp/4yysq1322155631.ps tmp/4yysq1322155631.png",intern=TRUE))
character(0)
> try(system("convert tmp/5otoq1322155631.ps tmp/5otoq1322155631.png",intern=TRUE))
character(0)
> try(system("convert tmp/66m0w1322155631.ps tmp/66m0w1322155631.png",intern=TRUE))
character(0)
> try(system("convert tmp/77j6x1322155631.ps tmp/77j6x1322155631.png",intern=TRUE))
character(0)
> try(system("convert tmp/8p3vh1322155631.ps tmp/8p3vh1322155631.png",intern=TRUE))
character(0)
> try(system("convert tmp/9s7zf1322155631.ps tmp/9s7zf1322155631.png",intern=TRUE))
character(0)
> try(system("convert tmp/10moey1322155631.ps tmp/10moey1322155631.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
2.790 0.260 3.082