R version 2.13.0 (2011-04-13)
Copyright (C) 2011 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(170588
+ ,95556
+ ,114468
+ ,86621
+ ,54565
+ ,88594
+ ,113337
+ ,63016
+ ,74151
+ ,152510
+ ,79774
+ ,77921
+ ,86206
+ ,31258
+ ,53212
+ ,37257
+ ,52491
+ ,34956
+ ,306055
+ ,91256
+ ,149703
+ ,32750
+ ,22807
+ ,6853
+ ,116502
+ ,77411
+ ,58907
+ ,130539
+ ,48821
+ ,67067
+ ,161876
+ ,52295
+ ,110563
+ ,128274
+ ,63262
+ ,58126
+ ,102350
+ ,50466
+ ,57113
+ ,193024
+ ,62932
+ ,77993
+ ,141574
+ ,38439
+ ,68091
+ ,253559
+ ,70817
+ ,124676
+ ,181110
+ ,105965
+ ,109522
+ ,198432
+ ,73795
+ ,75865
+ ,113853
+ ,82043
+ ,79746
+ ,159940
+ ,74349
+ ,77844
+ ,166822
+ ,82204
+ ,98681
+ ,286675
+ ,55709
+ ,105531
+ ,91657
+ ,37137
+ ,51428
+ ,108278
+ ,70780
+ ,65703
+ ,146342
+ ,55027
+ ,72562
+ ,145142
+ ,56699
+ ,81728
+ ,161740
+ ,65911
+ ,95580
+ ,160905
+ ,56316
+ ,98278
+ ,106888
+ ,26982
+ ,46629
+ ,188150
+ ,54628
+ ,115189
+ ,189401
+ ,96750
+ ,124865
+ ,129484
+ ,53009
+ ,59392
+ ,204030
+ ,64664
+ ,127818
+ ,68538
+ ,36990
+ ,17821
+ ,243625
+ ,85224
+ ,154076
+ ,167255
+ ,37048
+ ,64881
+ ,264528
+ ,59635
+ ,136506
+ ,122024
+ ,42051
+ ,66524
+ ,80964
+ ,26998
+ ,45988
+ ,209795
+ ,63717
+ ,107445
+ ,224205
+ ,55071
+ ,102772
+ ,115971
+ ,40001
+ ,46657
+ ,138191
+ ,54506
+ ,97563
+ ,81106
+ ,35838
+ ,36663
+ ,93125
+ ,50838
+ ,55369
+ ,305756
+ ,86997
+ ,77921
+ ,78800
+ ,33032
+ ,56968
+ ,158835
+ ,61704
+ ,77519
+ ,221745
+ ,117986
+ ,129805
+ ,131108
+ ,56733
+ ,72761
+ ,128734
+ ,55064
+ ,81278
+ ,24188
+ ,5950
+ ,15049
+ ,257662
+ ,84607
+ ,113935
+ ,65029
+ ,32551
+ ,25109
+ ,98066
+ ,31701
+ ,45824
+ ,173587
+ ,71170
+ ,89644
+ ,180042
+ ,101773
+ ,109011
+ ,197266
+ ,101653
+ ,134245
+ ,212060
+ ,81493
+ ,136692
+ ,141582
+ ,55901
+ ,50741)
+ ,dim=c(3
+ ,60)
+ ,dimnames=list(c('TotalTime'
+ ,'Compendiumwriting_char'
+ ,'Compendiumwriting_sec')
+ ,1:60))
> y <- array(NA,dim=c(3,60),dimnames=list(c('TotalTime','Compendiumwriting_char','Compendiumwriting_sec'),1:60))
> 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'
> 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
TotalTime Compendiumwriting_char Compendiumwriting_sec
1 170588 95556 114468
2 86621 54565 88594
3 113337 63016 74151
4 152510 79774 77921
5 86206 31258 53212
6 37257 52491 34956
7 306055 91256 149703
8 32750 22807 6853
9 116502 77411 58907
10 130539 48821 67067
11 161876 52295 110563
12 128274 63262 58126
13 102350 50466 57113
14 193024 62932 77993
15 141574 38439 68091
16 253559 70817 124676
17 181110 105965 109522
18 198432 73795 75865
19 113853 82043 79746
20 159940 74349 77844
21 166822 82204 98681
22 286675 55709 105531
23 91657 37137 51428
24 108278 70780 65703
25 146342 55027 72562
26 145142 56699 81728
27 161740 65911 95580
28 160905 56316 98278
29 106888 26982 46629
30 188150 54628 115189
31 189401 96750 124865
32 129484 53009 59392
33 204030 64664 127818
34 68538 36990 17821
35 243625 85224 154076
36 167255 37048 64881
37 264528 59635 136506
38 122024 42051 66524
39 80964 26998 45988
40 209795 63717 107445
41 224205 55071 102772
42 115971 40001 46657
43 138191 54506 97563
44 81106 35838 36663
45 93125 50838 55369
46 305756 86997 77921
47 78800 33032 56968
48 158835 61704 77519
49 221745 117986 129805
50 131108 56733 72761
51 128734 55064 81278
52 24188 5950 15049
53 257662 84607 113935
54 65029 32551 25109
55 98066 31701 45824
56 173587 71170 89644
57 180042 101773 109011
58 197266 101653 134245
59 212060 81493 136692
60 141582 55901 50741
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Compendiumwriting_char Compendiumwriting_sec
2.065e+04 2.195e-01 1.461e+00
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-75484 -21247 -7313 12369 152132
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.065e+04 1.368e+04 1.510 0.137
Compendiumwriting_char 2.195e-01 3.097e-01 0.709 0.481
Compendiumwriting_sec 1.462e+00 2.075e-01 7.045 2.68e-09 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 36490 on 57 degrees of freedom
Multiple R-squared: 0.6937, Adjusted R-squared: 0.683
F-statistic: 64.56 on 2 and 57 DF, p-value: 2.257e-15
> 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.46002015 0.9200403 0.53997985
[2,] 0.70524379 0.5895124 0.29475621
[3,] 0.73499771 0.5300046 0.26500229
[4,] 0.63704727 0.7259055 0.36295273
[5,] 0.55176142 0.8964772 0.44823858
[6,] 0.47041525 0.9408305 0.52958475
[7,] 0.40965565 0.8193113 0.59034435
[8,] 0.31646538 0.6329308 0.68353462
[9,] 0.44594574 0.8918915 0.55405426
[10,] 0.39932969 0.7986594 0.60067031
[11,] 0.39289699 0.7857940 0.60710301
[12,] 0.33733738 0.6746748 0.66266262
[13,] 0.44732114 0.8946423 0.55267886
[14,] 0.45490143 0.9098029 0.54509857
[15,] 0.38339677 0.7667935 0.61660323
[16,] 0.31928004 0.6385601 0.68071996
[17,] 0.76278897 0.4744221 0.23721103
[18,] 0.70357491 0.5928502 0.29642509
[19,] 0.66831827 0.6633635 0.33168173
[20,] 0.59632945 0.8073411 0.40367055
[21,] 0.52368385 0.9526323 0.47631615
[22,] 0.46063735 0.9212747 0.53936265
[23,] 0.40643632 0.8128726 0.59356368
[24,] 0.34207919 0.6841584 0.65792081
[25,] 0.29451594 0.5890319 0.70548406
[26,] 0.29995098 0.5999020 0.70004902
[27,] 0.24331484 0.4866297 0.75668516
[28,] 0.20214770 0.4042954 0.79785230
[29,] 0.16174089 0.3234818 0.83825911
[30,] 0.12838040 0.2567608 0.87161960
[31,] 0.14397734 0.2879547 0.85602266
[32,] 0.17327886 0.3465577 0.82672114
[33,] 0.12790728 0.2558146 0.87209272
[34,] 0.09450323 0.1890065 0.90549677
[35,] 0.08233113 0.1646623 0.91766887
[36,] 0.17337345 0.3467469 0.82662655
[37,] 0.13282941 0.2656588 0.86717059
[38,] 0.11228485 0.2245697 0.88771515
[39,] 0.07857824 0.1571565 0.92142176
[40,] 0.06774183 0.1354837 0.93225817
[41,] 0.91629609 0.1674078 0.08370391
[42,] 0.89779328 0.2044134 0.10220672
[43,] 0.85114041 0.2977192 0.14885959
[44,] 0.78786934 0.4242613 0.21213066
[45,] 0.69180006 0.6163999 0.30819994
[46,] 0.59684432 0.8063114 0.40315568
[47,] 0.54223305 0.9155339 0.45776695
[48,] 0.89630700 0.2073860 0.10369300
[49,] 0.88897948 0.2220410 0.11102052
> postscript(file="/var/wessaorg/rcomp/tmp/1exi31321895392.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/wessaorg/rcomp/tmp/2rz2v1321895392.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/wessaorg/rcomp/tmp/36pko1321895392.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/wessaorg/rcomp/tmp/4w16m1321895392.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/wessaorg/rcomp/tmp/5fspd1321895392.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 = 60
Frequency = 1
1 2 3 4 5 6
-38327.2718 -75483.9239 -29514.5472 470.8767 -19073.9107 -46002.1463
7 8 9 10 11 12
46588.3556 -2922.6748 -7230.0359 1155.8658 -31837.8834 8788.6036
13 14 15 16 17 18
-12846.6324 44575.9004 12972.8135 35154.3644 -22861.2395 50709.8581
19 20 21 22 23 24
-41351.2843 9204.0145 -16090.6506 99566.0019 -12305.8809 -23930.9407
25 26 27 28 29 30
7566.0525 -7396.7764 -13064.8531 -15737.1347 12167.4178 -12836.6822
31 32 33 34 35 36
-34971.2717 10398.5631 -17616.2149 13723.1763 -20908.8551 43650.4262
37 38 39 40 41 42
31288.1845 -5079.7669 -12823.2881 18131.2530 41268.2313 18352.2658
43 44 45 46 47 48
-37008.9461 -993.1136 -19604.4611 152131.6722 -32358.5445 11349.1442
49 50 51 52 53 54
-14507.5435 -8333.1905 -22788.2853 -19763.3832 51928.6300 537.1559
55 56 57 58 59 60
3486.2442 6303.2907 -22262.4233 -41890.9516 -26248.7424 34505.0885
> postscript(file="/var/wessaorg/rcomp/tmp/6dlkv1321895392.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 = 60
Frequency = 1
lag(myerror, k = 1) myerror
0 -38327.2718 NA
1 -75483.9239 -38327.2718
2 -29514.5472 -75483.9239
3 470.8767 -29514.5472
4 -19073.9107 470.8767
5 -46002.1463 -19073.9107
6 46588.3556 -46002.1463
7 -2922.6748 46588.3556
8 -7230.0359 -2922.6748
9 1155.8658 -7230.0359
10 -31837.8834 1155.8658
11 8788.6036 -31837.8834
12 -12846.6324 8788.6036
13 44575.9004 -12846.6324
14 12972.8135 44575.9004
15 35154.3644 12972.8135
16 -22861.2395 35154.3644
17 50709.8581 -22861.2395
18 -41351.2843 50709.8581
19 9204.0145 -41351.2843
20 -16090.6506 9204.0145
21 99566.0019 -16090.6506
22 -12305.8809 99566.0019
23 -23930.9407 -12305.8809
24 7566.0525 -23930.9407
25 -7396.7764 7566.0525
26 -13064.8531 -7396.7764
27 -15737.1347 -13064.8531
28 12167.4178 -15737.1347
29 -12836.6822 12167.4178
30 -34971.2717 -12836.6822
31 10398.5631 -34971.2717
32 -17616.2149 10398.5631
33 13723.1763 -17616.2149
34 -20908.8551 13723.1763
35 43650.4262 -20908.8551
36 31288.1845 43650.4262
37 -5079.7669 31288.1845
38 -12823.2881 -5079.7669
39 18131.2530 -12823.2881
40 41268.2313 18131.2530
41 18352.2658 41268.2313
42 -37008.9461 18352.2658
43 -993.1136 -37008.9461
44 -19604.4611 -993.1136
45 152131.6722 -19604.4611
46 -32358.5445 152131.6722
47 11349.1442 -32358.5445
48 -14507.5435 11349.1442
49 -8333.1905 -14507.5435
50 -22788.2853 -8333.1905
51 -19763.3832 -22788.2853
52 51928.6300 -19763.3832
53 537.1559 51928.6300
54 3486.2442 537.1559
55 6303.2907 3486.2442
56 -22262.4233 6303.2907
57 -41890.9516 -22262.4233
58 -26248.7424 -41890.9516
59 34505.0885 -26248.7424
60 NA 34505.0885
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -75483.9239 -38327.2718
[2,] -29514.5472 -75483.9239
[3,] 470.8767 -29514.5472
[4,] -19073.9107 470.8767
[5,] -46002.1463 -19073.9107
[6,] 46588.3556 -46002.1463
[7,] -2922.6748 46588.3556
[8,] -7230.0359 -2922.6748
[9,] 1155.8658 -7230.0359
[10,] -31837.8834 1155.8658
[11,] 8788.6036 -31837.8834
[12,] -12846.6324 8788.6036
[13,] 44575.9004 -12846.6324
[14,] 12972.8135 44575.9004
[15,] 35154.3644 12972.8135
[16,] -22861.2395 35154.3644
[17,] 50709.8581 -22861.2395
[18,] -41351.2843 50709.8581
[19,] 9204.0145 -41351.2843
[20,] -16090.6506 9204.0145
[21,] 99566.0019 -16090.6506
[22,] -12305.8809 99566.0019
[23,] -23930.9407 -12305.8809
[24,] 7566.0525 -23930.9407
[25,] -7396.7764 7566.0525
[26,] -13064.8531 -7396.7764
[27,] -15737.1347 -13064.8531
[28,] 12167.4178 -15737.1347
[29,] -12836.6822 12167.4178
[30,] -34971.2717 -12836.6822
[31,] 10398.5631 -34971.2717
[32,] -17616.2149 10398.5631
[33,] 13723.1763 -17616.2149
[34,] -20908.8551 13723.1763
[35,] 43650.4262 -20908.8551
[36,] 31288.1845 43650.4262
[37,] -5079.7669 31288.1845
[38,] -12823.2881 -5079.7669
[39,] 18131.2530 -12823.2881
[40,] 41268.2313 18131.2530
[41,] 18352.2658 41268.2313
[42,] -37008.9461 18352.2658
[43,] -993.1136 -37008.9461
[44,] -19604.4611 -993.1136
[45,] 152131.6722 -19604.4611
[46,] -32358.5445 152131.6722
[47,] 11349.1442 -32358.5445
[48,] -14507.5435 11349.1442
[49,] -8333.1905 -14507.5435
[50,] -22788.2853 -8333.1905
[51,] -19763.3832 -22788.2853
[52,] 51928.6300 -19763.3832
[53,] 537.1559 51928.6300
[54,] 3486.2442 537.1559
[55,] 6303.2907 3486.2442
[56,] -22262.4233 6303.2907
[57,] -41890.9516 -22262.4233
[58,] -26248.7424 -41890.9516
[59,] 34505.0885 -26248.7424
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -75483.9239 -38327.2718
2 -29514.5472 -75483.9239
3 470.8767 -29514.5472
4 -19073.9107 470.8767
5 -46002.1463 -19073.9107
6 46588.3556 -46002.1463
7 -2922.6748 46588.3556
8 -7230.0359 -2922.6748
9 1155.8658 -7230.0359
10 -31837.8834 1155.8658
11 8788.6036 -31837.8834
12 -12846.6324 8788.6036
13 44575.9004 -12846.6324
14 12972.8135 44575.9004
15 35154.3644 12972.8135
16 -22861.2395 35154.3644
17 50709.8581 -22861.2395
18 -41351.2843 50709.8581
19 9204.0145 -41351.2843
20 -16090.6506 9204.0145
21 99566.0019 -16090.6506
22 -12305.8809 99566.0019
23 -23930.9407 -12305.8809
24 7566.0525 -23930.9407
25 -7396.7764 7566.0525
26 -13064.8531 -7396.7764
27 -15737.1347 -13064.8531
28 12167.4178 -15737.1347
29 -12836.6822 12167.4178
30 -34971.2717 -12836.6822
31 10398.5631 -34971.2717
32 -17616.2149 10398.5631
33 13723.1763 -17616.2149
34 -20908.8551 13723.1763
35 43650.4262 -20908.8551
36 31288.1845 43650.4262
37 -5079.7669 31288.1845
38 -12823.2881 -5079.7669
39 18131.2530 -12823.2881
40 41268.2313 18131.2530
41 18352.2658 41268.2313
42 -37008.9461 18352.2658
43 -993.1136 -37008.9461
44 -19604.4611 -993.1136
45 152131.6722 -19604.4611
46 -32358.5445 152131.6722
47 11349.1442 -32358.5445
48 -14507.5435 11349.1442
49 -8333.1905 -14507.5435
50 -22788.2853 -8333.1905
51 -19763.3832 -22788.2853
52 51928.6300 -19763.3832
53 537.1559 51928.6300
54 3486.2442 537.1559
55 6303.2907 3486.2442
56 -22262.4233 6303.2907
57 -41890.9516 -22262.4233
58 -26248.7424 -41890.9516
59 34505.0885 -26248.7424
> 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/wessaorg/rcomp/tmp/7z2hl1321895392.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/wessaorg/rcomp/tmp/8deps1321895392.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/wessaorg/rcomp/tmp/9c8w31321895392.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/wessaorg/rcomp/tmp/107yyd1321895392.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/wessaorg/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/wessaorg/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/wessaorg/rcomp/tmp/11ni2j1321895392.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/wessaorg/rcomp/tmp/12k86z1321895392.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/wessaorg/rcomp/tmp/131wcf1321895392.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/wessaorg/rcomp/tmp/14ixah1321895392.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/wessaorg/rcomp/tmp/15zs3d1321895392.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/wessaorg/rcomp/tmp/166s5j1321895392.tab")
+ }
>
> try(system("convert tmp/1exi31321895392.ps tmp/1exi31321895392.png",intern=TRUE))
character(0)
> try(system("convert tmp/2rz2v1321895392.ps tmp/2rz2v1321895392.png",intern=TRUE))
character(0)
> try(system("convert tmp/36pko1321895392.ps tmp/36pko1321895392.png",intern=TRUE))
character(0)
> try(system("convert tmp/4w16m1321895392.ps tmp/4w16m1321895392.png",intern=TRUE))
character(0)
> try(system("convert tmp/5fspd1321895392.ps tmp/5fspd1321895392.png",intern=TRUE))
character(0)
> try(system("convert tmp/6dlkv1321895392.ps tmp/6dlkv1321895392.png",intern=TRUE))
character(0)
> try(system("convert tmp/7z2hl1321895392.ps tmp/7z2hl1321895392.png",intern=TRUE))
character(0)
> try(system("convert tmp/8deps1321895392.ps tmp/8deps1321895392.png",intern=TRUE))
character(0)
> try(system("convert tmp/9c8w31321895392.ps tmp/9c8w31321895392.png",intern=TRUE))
character(0)
> try(system("convert tmp/107yyd1321895392.ps tmp/107yyd1321895392.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.394 0.476 3.935