R version 2.15.2 (2012-10-26) -- "Trick or Treat"
Copyright (C) 2012 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: i686-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(8500
+ ,7300
+ ,8000
+ ,8350
+ ,7300
+ ,7700
+ ,8300
+ ,7200
+ ,9000
+ ,8400
+ ,8500
+ ,9200
+ ,9000
+ ,10000
+ ,9500
+ ,8300
+ ,9200
+ ,9200
+ ,7000
+ ,8200
+ ,8900
+ ,10300
+ ,8000
+ ,8700
+ ,7150
+ ,9000
+ ,10000
+ ,8100
+ ,9500
+ ,10500
+ ,7200
+ ,7400
+ ,9800
+ ,6000
+ ,10500
+ ,9900
+ ,6750
+ ,7600
+ ,10000
+ ,9200
+ ,8300
+ ,8900
+ ,7600
+ ,8400
+ ,10000
+ ,7000
+ ,4300
+ ,6500
+ ,8288
+ ,7400
+ ,9000
+ ,8400
+ ,8200
+ ,8700
+ ,14000
+ ,7500
+ ,9500
+ ,8500
+ ,8500
+ ,10000
+ ,9500
+ ,7600
+ ,8900
+ ,11811
+ ,8200
+ ,8700
+ ,10000
+ ,9200
+ ,9800
+ ,9500
+ ,7500
+ ,8750
+ ,9500
+ ,8200
+ ,9000
+ ,9452
+ ,9500
+ ,7500
+ ,9500
+ ,10000
+ ,7900
+ ,8600
+ ,11000
+ ,8700
+ ,11763
+ ,7500
+ ,8800
+ ,9766
+ ,9000
+ ,8900
+ ,11400
+ ,8000
+ ,9500
+ ,9500
+ ,8900
+ ,9600
+ ,11994
+ ,8700
+ ,8700
+ ,8400
+ ,6000
+ ,10000
+ ,7360
+ ,5000
+ ,10500
+ ,7400
+ ,7500
+ ,9250
+ ,8558
+ ,8000
+ ,8250
+ ,7000
+ ,7600
+ ,6500
+ ,7400
+ ,8900
+ ,8525
+ ,7200
+ ,9200
+ ,9255
+ ,8600
+ ,9500
+ ,8750
+ ,7800
+ ,7800
+ ,8540
+ ,7500
+ ,10000
+ ,7890
+ ,9000
+ ,8500
+ ,9000
+ ,7429
+ ,9000
+ ,10250
+ ,7206
+ ,9200
+ ,10100
+ ,7613
+ ,9450
+ ,8560
+ ,7200
+ ,8600
+ ,8000
+ ,7500
+ ,10500
+ ,7800
+ ,7500
+ ,7700
+ ,8000
+ ,9071
+ ,8500
+ ,1000
+ ,7600
+ ,7000
+ ,10250
+ ,8359
+ ,7500
+ ,9500
+ ,15000
+ ,6500
+ ,9350
+ ,6500
+ ,9000
+ ,8900
+ ,6500
+ ,8500
+ ,8750
+ ,6125
+ ,8600
+ ,8250
+ ,6000
+ ,9400
+ ,7900
+ ,7000
+ ,9200
+ ,9000
+ ,8500
+ ,8400
+ ,9250
+ ,7000
+ ,7900
+ ,10000
+ ,7000
+ ,10000
+ ,10250
+ ,6600
+ ,9000
+ ,9500
+ ,6800
+ ,8500
+ ,8900
+ ,12000
+ ,11000
+ ,9750
+ ,7200
+ ,9500
+ ,9000
+ ,7200
+ ,8900
+ ,8750
+ ,7300
+ ,9000
+ ,8000
+ ,7500
+ ,9200
+ ,8450
+ ,7000
+ ,8700
+ ,8700
+ ,7000
+ ,7000
+ ,8500
+ ,6000
+ ,7500
+ ,7500)
+ ,dim=c(3
+ ,72)
+ ,dimnames=list(c('A'
+ ,'B'
+ ,'C')
+ ,1:72))
> y <- array(NA,dim=c(3,72),dimnames=list(c('A','B','C'),1:72))
> 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 = 'Do not include Seasonal Dummies'
> par1 = '1'
> par3 <- '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
Attaching package: 'zoo'
The following object(s) are masked from 'package:base':
as.Date, 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
A B C t
1 8500 7300 8000 1
2 8350 7300 7700 2
3 8300 7200 9000 3
4 8400 8500 9200 4
5 9000 10000 9500 5
6 8300 9200 9200 6
7 7000 8200 8900 7
8 10300 8000 8700 8
9 7150 9000 10000 9
10 8100 9500 10500 10
11 7200 7400 9800 11
12 6000 10500 9900 12
13 6750 7600 10000 13
14 9200 8300 8900 14
15 7600 8400 10000 15
16 7000 4300 6500 16
17 8288 7400 9000 17
18 8400 8200 8700 18
19 14000 7500 9500 19
20 8500 8500 10000 20
21 9500 7600 8900 21
22 11811 8200 8700 22
23 10000 9200 9800 23
24 9500 7500 8750 24
25 9500 8200 9000 25
26 9452 9500 7500 26
27 9500 10000 7900 27
28 8600 11000 8700 28
29 11763 7500 8800 29
30 9766 9000 8900 30
31 11400 8000 9500 31
32 9500 8900 9600 32
33 11994 8700 8700 33
34 8400 6000 10000 34
35 7360 5000 10500 35
36 7400 7500 9250 36
37 8558 8000 8250 37
38 7000 7600 6500 38
39 7400 8900 8525 39
40 7200 9200 9255 40
41 8600 9500 8750 41
42 7800 7800 8540 42
43 7500 10000 7890 43
44 9000 8500 9000 44
45 7429 9000 10250 45
46 7206 9200 10100 46
47 7613 9450 8560 47
48 7200 8600 8000 48
49 7500 10500 7800 49
50 7500 7700 8000 50
51 9071 8500 1000 51
52 7600 7000 10250 52
53 8359 7500 9500 53
54 15000 6500 9350 54
55 6500 9000 8900 55
56 6500 8500 8750 56
57 6125 8600 8250 57
58 6000 9400 7900 58
59 7000 9200 9000 59
60 8500 8400 9250 60
61 7000 7900 10000 61
62 7000 10000 10250 62
63 6600 9000 9500 63
64 6800 8500 8900 64
65 12000 11000 9750 65
66 7200 9500 9000 66
67 7200 8900 8750 67
68 7300 9000 8000 68
69 7500 9200 8450 69
70 7000 8700 8700 70
71 7000 7000 8500 71
72 6000 7500 7500 72
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) B C t
9.530e+03 -5.910e-02 5.226e-03 -2.095e+01
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-2709.4 -885.2 -487.7 709.7 6936.8
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 9.530e+03 2.050e+03 4.649 1.59e-05 ***
B -5.910e-02 1.762e-01 -0.335 0.7384
C 5.226e-03 1.661e-01 0.031 0.9750
t -2.095e+01 1.020e+01 -2.054 0.0438 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1755 on 68 degrees of freedom
Multiple R-squared: 0.06599, Adjusted R-squared: 0.02478
F-statistic: 1.601 on 3 and 68 DF, p-value: 0.1971
> 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.002161744 0.004323489 0.99783826
[2,] 0.196805466 0.393610931 0.80319453
[3,] 0.146116727 0.292233455 0.85388327
[4,] 0.079608009 0.159216017 0.92039199
[5,] 0.043792765 0.087585529 0.95620724
[6,] 0.083692764 0.167385528 0.91630724
[7,] 0.055843865 0.111687730 0.94415613
[8,] 0.059096343 0.118192686 0.94090366
[9,] 0.037966098 0.075932195 0.96203390
[10,] 0.039476609 0.078953219 0.96052339
[11,] 0.031692307 0.063384615 0.96830769
[12,] 0.021923672 0.043847344 0.97807633
[13,] 0.727259305 0.545481391 0.27274070
[14,] 0.668062515 0.663874970 0.33193749
[15,] 0.598428047 0.803143906 0.40157195
[16,] 0.652815588 0.694368825 0.34718441
[17,] 0.581108091 0.837783818 0.41889191
[18,] 0.506318905 0.987362190 0.49368109
[19,] 0.433224481 0.866448963 0.56677552
[20,] 0.377243227 0.754486453 0.62275677
[21,] 0.313349987 0.626699974 0.68665001
[22,] 0.270322569 0.540645139 0.72967743
[23,] 0.302644726 0.605289451 0.69735527
[24,] 0.248139723 0.496279446 0.75186028
[25,] 0.259265245 0.518530491 0.74073475
[26,] 0.220117777 0.440235555 0.77988222
[27,] 0.323135146 0.646270292 0.67686485
[28,] 0.329584615 0.659169230 0.67041538
[29,] 0.348497320 0.696994640 0.65150268
[30,] 0.377344160 0.754688320 0.62265584
[31,] 0.347562225 0.695124449 0.65243778
[32,] 0.415349950 0.830699899 0.58465005
[33,] 0.408355864 0.816711728 0.59164414
[34,] 0.400601165 0.801202329 0.59939884
[35,] 0.341947800 0.683895599 0.65805220
[36,] 0.297565926 0.595131853 0.70243407
[37,] 0.260383276 0.520766552 0.73961672
[38,] 0.210606643 0.421213285 0.78939336
[39,] 0.177836843 0.355673686 0.82216316
[40,] 0.151650717 0.303301433 0.84834928
[41,] 0.118808979 0.237617958 0.88119102
[42,] 0.099703876 0.199407752 0.90029612
[43,] 0.074523643 0.149047286 0.92547636
[44,] 0.057610200 0.115220400 0.94238980
[45,] 0.069871065 0.139742130 0.93012893
[46,] 0.069627588 0.139255176 0.93037241
[47,] 0.049834659 0.099669317 0.95016534
[48,] 0.983331427 0.033337146 0.01666857
[49,] 0.973362514 0.053274972 0.02663749
[50,] 0.958973338 0.082053323 0.04102666
[51,] 0.937633357 0.124733287 0.06236664
[52,] 0.911230302 0.177539395 0.08876970
[53,] 0.864924871 0.270150258 0.13507513
[54,] 0.908020128 0.183959745 0.09197987
[55,] 0.899615756 0.200768487 0.10038424
[56,] 0.918458868 0.163082265 0.08154113
[57,] 0.883427434 0.233145133 0.11657257
[58,] 0.784968275 0.430063451 0.21503173
[59,] 0.981728474 0.036543052 0.01827153
> postscript(file="/var/wessaorg/rcomp/tmp/1r9jb1352122008.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/2n24d1352122008.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/3ys6t1352122008.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/4wsr61352122008.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/5oo5g1352122008.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 = 72
Frequency = 1
1 2 3 4 5 6 7
-618.9915 -746.4783 -788.2371 -591.5088 116.5165 -608.2492 -1944.8347
8 9 10 11 12 13 14
1365.3361 -1711.4142 -713.5328 -1713.0362 -2709.4078 -2110.3712 407.6921
15 16 17 18 19 20 21
-1171.2017 -1974.2687 -495.1835 -313.3914 5262.0038 -160.5654 812.9400
22 23 24 25 26 27 28
3181.3898 1444.6847 870.6499 931.6577 989.2707 1085.6748 261.5376
29 30 31 32 33 34 35
3238.1150 1350.1856 2942.8965 1116.5079 3624.3371 -115.0781 -1195.8447
36 37 38 39 40 41 42
-980.6199 233.1011 -1318.4470 -831.2568 -996.3972 444.9170 -433.5079
43 44 45 46 47 48 49
-579.1485 847.3476 -679.6907 -869.1418 -418.3733 -857.7351 -423.4571
50 51 52 53 54 55 56
-569.0333 1106.7752 -480.2710 333.1433 6936.7739 -1392.1824 -1400.0025
57 58 59 60 61 62 63
-1745.5342 -1800.4808 -797.1042 675.2556 -837.2682 -693.5223 -1127.7559
64 65 66 67 68 69 70
-933.2242 4431.0253 -432.7576 -445.9649 -315.1901 -84.7769 -594.6875
71 72
-673.1647 -1617.4437
> postscript(file="/var/wessaorg/rcomp/tmp/6ebur1352122008.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 = 72
Frequency = 1
lag(myerror, k = 1) myerror
0 -618.9915 NA
1 -746.4783 -618.9915
2 -788.2371 -746.4783
3 -591.5088 -788.2371
4 116.5165 -591.5088
5 -608.2492 116.5165
6 -1944.8347 -608.2492
7 1365.3361 -1944.8347
8 -1711.4142 1365.3361
9 -713.5328 -1711.4142
10 -1713.0362 -713.5328
11 -2709.4078 -1713.0362
12 -2110.3712 -2709.4078
13 407.6921 -2110.3712
14 -1171.2017 407.6921
15 -1974.2687 -1171.2017
16 -495.1835 -1974.2687
17 -313.3914 -495.1835
18 5262.0038 -313.3914
19 -160.5654 5262.0038
20 812.9400 -160.5654
21 3181.3898 812.9400
22 1444.6847 3181.3898
23 870.6499 1444.6847
24 931.6577 870.6499
25 989.2707 931.6577
26 1085.6748 989.2707
27 261.5376 1085.6748
28 3238.1150 261.5376
29 1350.1856 3238.1150
30 2942.8965 1350.1856
31 1116.5079 2942.8965
32 3624.3371 1116.5079
33 -115.0781 3624.3371
34 -1195.8447 -115.0781
35 -980.6199 -1195.8447
36 233.1011 -980.6199
37 -1318.4470 233.1011
38 -831.2568 -1318.4470
39 -996.3972 -831.2568
40 444.9170 -996.3972
41 -433.5079 444.9170
42 -579.1485 -433.5079
43 847.3476 -579.1485
44 -679.6907 847.3476
45 -869.1418 -679.6907
46 -418.3733 -869.1418
47 -857.7351 -418.3733
48 -423.4571 -857.7351
49 -569.0333 -423.4571
50 1106.7752 -569.0333
51 -480.2710 1106.7752
52 333.1433 -480.2710
53 6936.7739 333.1433
54 -1392.1824 6936.7739
55 -1400.0025 -1392.1824
56 -1745.5342 -1400.0025
57 -1800.4808 -1745.5342
58 -797.1042 -1800.4808
59 675.2556 -797.1042
60 -837.2682 675.2556
61 -693.5223 -837.2682
62 -1127.7559 -693.5223
63 -933.2242 -1127.7559
64 4431.0253 -933.2242
65 -432.7576 4431.0253
66 -445.9649 -432.7576
67 -315.1901 -445.9649
68 -84.7769 -315.1901
69 -594.6875 -84.7769
70 -673.1647 -594.6875
71 -1617.4437 -673.1647
72 NA -1617.4437
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -746.4783 -618.9915
[2,] -788.2371 -746.4783
[3,] -591.5088 -788.2371
[4,] 116.5165 -591.5088
[5,] -608.2492 116.5165
[6,] -1944.8347 -608.2492
[7,] 1365.3361 -1944.8347
[8,] -1711.4142 1365.3361
[9,] -713.5328 -1711.4142
[10,] -1713.0362 -713.5328
[11,] -2709.4078 -1713.0362
[12,] -2110.3712 -2709.4078
[13,] 407.6921 -2110.3712
[14,] -1171.2017 407.6921
[15,] -1974.2687 -1171.2017
[16,] -495.1835 -1974.2687
[17,] -313.3914 -495.1835
[18,] 5262.0038 -313.3914
[19,] -160.5654 5262.0038
[20,] 812.9400 -160.5654
[21,] 3181.3898 812.9400
[22,] 1444.6847 3181.3898
[23,] 870.6499 1444.6847
[24,] 931.6577 870.6499
[25,] 989.2707 931.6577
[26,] 1085.6748 989.2707
[27,] 261.5376 1085.6748
[28,] 3238.1150 261.5376
[29,] 1350.1856 3238.1150
[30,] 2942.8965 1350.1856
[31,] 1116.5079 2942.8965
[32,] 3624.3371 1116.5079
[33,] -115.0781 3624.3371
[34,] -1195.8447 -115.0781
[35,] -980.6199 -1195.8447
[36,] 233.1011 -980.6199
[37,] -1318.4470 233.1011
[38,] -831.2568 -1318.4470
[39,] -996.3972 -831.2568
[40,] 444.9170 -996.3972
[41,] -433.5079 444.9170
[42,] -579.1485 -433.5079
[43,] 847.3476 -579.1485
[44,] -679.6907 847.3476
[45,] -869.1418 -679.6907
[46,] -418.3733 -869.1418
[47,] -857.7351 -418.3733
[48,] -423.4571 -857.7351
[49,] -569.0333 -423.4571
[50,] 1106.7752 -569.0333
[51,] -480.2710 1106.7752
[52,] 333.1433 -480.2710
[53,] 6936.7739 333.1433
[54,] -1392.1824 6936.7739
[55,] -1400.0025 -1392.1824
[56,] -1745.5342 -1400.0025
[57,] -1800.4808 -1745.5342
[58,] -797.1042 -1800.4808
[59,] 675.2556 -797.1042
[60,] -837.2682 675.2556
[61,] -693.5223 -837.2682
[62,] -1127.7559 -693.5223
[63,] -933.2242 -1127.7559
[64,] 4431.0253 -933.2242
[65,] -432.7576 4431.0253
[66,] -445.9649 -432.7576
[67,] -315.1901 -445.9649
[68,] -84.7769 -315.1901
[69,] -594.6875 -84.7769
[70,] -673.1647 -594.6875
[71,] -1617.4437 -673.1647
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -746.4783 -618.9915
2 -788.2371 -746.4783
3 -591.5088 -788.2371
4 116.5165 -591.5088
5 -608.2492 116.5165
6 -1944.8347 -608.2492
7 1365.3361 -1944.8347
8 -1711.4142 1365.3361
9 -713.5328 -1711.4142
10 -1713.0362 -713.5328
11 -2709.4078 -1713.0362
12 -2110.3712 -2709.4078
13 407.6921 -2110.3712
14 -1171.2017 407.6921
15 -1974.2687 -1171.2017
16 -495.1835 -1974.2687
17 -313.3914 -495.1835
18 5262.0038 -313.3914
19 -160.5654 5262.0038
20 812.9400 -160.5654
21 3181.3898 812.9400
22 1444.6847 3181.3898
23 870.6499 1444.6847
24 931.6577 870.6499
25 989.2707 931.6577
26 1085.6748 989.2707
27 261.5376 1085.6748
28 3238.1150 261.5376
29 1350.1856 3238.1150
30 2942.8965 1350.1856
31 1116.5079 2942.8965
32 3624.3371 1116.5079
33 -115.0781 3624.3371
34 -1195.8447 -115.0781
35 -980.6199 -1195.8447
36 233.1011 -980.6199
37 -1318.4470 233.1011
38 -831.2568 -1318.4470
39 -996.3972 -831.2568
40 444.9170 -996.3972
41 -433.5079 444.9170
42 -579.1485 -433.5079
43 847.3476 -579.1485
44 -679.6907 847.3476
45 -869.1418 -679.6907
46 -418.3733 -869.1418
47 -857.7351 -418.3733
48 -423.4571 -857.7351
49 -569.0333 -423.4571
50 1106.7752 -569.0333
51 -480.2710 1106.7752
52 333.1433 -480.2710
53 6936.7739 333.1433
54 -1392.1824 6936.7739
55 -1400.0025 -1392.1824
56 -1745.5342 -1400.0025
57 -1800.4808 -1745.5342
58 -797.1042 -1800.4808
59 675.2556 -797.1042
60 -837.2682 675.2556
61 -693.5223 -837.2682
62 -1127.7559 -693.5223
63 -933.2242 -1127.7559
64 4431.0253 -933.2242
65 -432.7576 4431.0253
66 -445.9649 -432.7576
67 -315.1901 -445.9649
68 -84.7769 -315.1901
69 -594.6875 -84.7769
70 -673.1647 -594.6875
71 -1617.4437 -673.1647
> 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/7ltjh1352122008.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/8ipge1352122008.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/94gqe1352122008.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/10pzw21352122008.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/11uwq01352122008.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/12128s1352122008.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/13ed3c1352122009.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/14b8g21352122009.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/150n5k1352122009.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/16t5l51352122009.tab")
+ }
>
> try(system("convert tmp/1r9jb1352122008.ps tmp/1r9jb1352122008.png",intern=TRUE))
character(0)
> try(system("convert tmp/2n24d1352122008.ps tmp/2n24d1352122008.png",intern=TRUE))
character(0)
> try(system("convert tmp/3ys6t1352122008.ps tmp/3ys6t1352122008.png",intern=TRUE))
character(0)
> try(system("convert tmp/4wsr61352122008.ps tmp/4wsr61352122008.png",intern=TRUE))
character(0)
> try(system("convert tmp/5oo5g1352122008.ps tmp/5oo5g1352122008.png",intern=TRUE))
character(0)
> try(system("convert tmp/6ebur1352122008.ps tmp/6ebur1352122008.png",intern=TRUE))
character(0)
> try(system("convert tmp/7ltjh1352122008.ps tmp/7ltjh1352122008.png",intern=TRUE))
character(0)
> try(system("convert tmp/8ipge1352122008.ps tmp/8ipge1352122008.png",intern=TRUE))
character(0)
> try(system("convert tmp/94gqe1352122008.ps tmp/94gqe1352122008.png",intern=TRUE))
character(0)
> try(system("convert tmp/10pzw21352122008.ps tmp/10pzw21352122008.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
6.560 0.966 7.522