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(458.15
+ ,000
+ ,477.59
+ ,050
+ ,504.91
+ ,100
+ ,502.61
+ ,150
+ ,514.12
+ ,200
+ ,512.64
+ ,250
+ ,546.06
+ ,300
+ ,525.26
+ ,350
+ ,571.20
+ ,400
+ ,551.22
+ ,450
+ ,604.26
+ ,500
+ ,510.28
+ ,550
+ ,574.91
+ ,600
+ ,580.80
+ ,650
+ ,527.33
+ ,700
+ ,571.37
+ ,750
+ ,587.97
+ ,800
+ ,557.65
+ ,850
+ ,619.61
+ ,900
+ ,631.11
+ ,950
+ ,583.14
+ ,1000
+ ,589.40
+ ,1050
+ ,603.19
+ ,1100
+ ,642.68
+ ,1150
+ ,615.91
+ ,1200
+ ,650.56
+ ,1250
+ ,607.41
+ ,1300
+ ,673.46
+ ,1350
+ ,680.11
+ ,1400
+ ,665.89
+ ,1450
+ ,711.79
+ ,1500
+ ,636.29
+ ,1550
+ ,580.08
+ ,1600
+ ,595.64
+ ,1650
+ ,661.80
+ ,1700
+ ,657.74
+ ,1750
+ ,646.05
+ ,1800
+ ,706.03
+ ,1850
+ ,712.38
+ ,1900
+ ,718.78
+ ,1950
+ ,699.49
+ ,2000
+ ,635.36
+ ,2050
+ ,682.09
+ ,2100
+ ,722.70
+ ,2150
+ ,731.22
+ ,2200
+ ,763.95
+ ,2250
+ ,739.86
+ ,2300
+ ,744.88
+ ,2350
+ ,746.73
+ ,2400
+ ,821.77
+ ,2450
+ ,752.76
+ ,2500
+ ,733.80
+ ,2550
+ ,735.91
+ ,2600
+ ,783.64
+ ,2650
+ ,711.28
+ ,2700
+ ,764.41
+ ,2750
+ ,833.71
+ ,2800
+ ,827.09
+ ,2850
+ ,766.46
+ ,2900
+ ,748.42
+ ,2950
+ ,870.61
+ ,3000
+ ,854.52
+ ,3050
+ ,858.34
+ ,3100
+ ,787.99
+ ,3150
+ ,834.26
+ ,3200
+ ,827.86
+ ,3250
+ ,771.05
+ ,3300
+ ,806.20
+ ,3350
+ ,873.40
+ ,3400
+ ,792.56
+ ,3450
+ ,855.02
+ ,3500
+ ,794.63
+ ,3550
+ ,861.60
+ ,3600
+ ,859.60
+ ,3650
+ ,856.97
+ ,3700
+ ,905.18
+ ,3750
+ ,933.00
+ ,3800
+ ,838.89
+ ,3850
+ ,903.42
+ ,3900
+ ,889.50
+ ,3950
+ ,914.18
+ ,4000
+ ,863.96
+ ,4050
+ ,937.39
+ ,4000
+ ,948.76
+ ,4150
+ ,900.66
+ ,4200
+ ,947.49
+ ,4250
+ ,904.22
+ ,4300
+ ,861.64
+ ,4350
+ ,918.50
+ ,4400
+ ,906.68
+ ,4450
+ ,966.54
+ ,4500
+ ,997.92
+ ,4550
+ ,965.90
+ ,4600
+ ,969.30
+ ,4650
+ ,904.80
+ ,4700
+ ,957.80
+ ,4750
+ ,1026.98
+ ,4800
+ ,1048.42
+ ,4850
+ ,953.19
+ ,4900
+ ,1020.25
+ ,4950)
+ ,dim=c(2
+ ,100)
+ ,dimnames=list(c('Oogst'
+ ,'Mest')
+ ,1:100))
> y <- array(NA,dim=c(2,100),dimnames=list(c('Oogst','Mest'),1:100))
> 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
Oogst Mest t
1 458.15 0 1
2 477.59 50 2
3 504.91 100 3
4 502.61 150 4
5 514.12 200 5
6 512.64 250 6
7 546.06 300 7
8 525.26 350 8
9 571.20 400 9
10 551.22 450 10
11 604.26 500 11
12 510.28 550 12
13 574.91 600 13
14 580.80 650 14
15 527.33 700 15
16 571.37 750 16
17 587.97 800 17
18 557.65 850 18
19 619.61 900 19
20 631.11 950 20
21 583.14 1000 21
22 589.40 1050 22
23 603.19 1100 23
24 642.68 1150 24
25 615.91 1200 25
26 650.56 1250 26
27 607.41 1300 27
28 673.46 1350 28
29 680.11 1400 29
30 665.89 1450 30
31 711.79 1500 31
32 636.29 1550 32
33 580.08 1600 33
34 595.64 1650 34
35 661.80 1700 35
36 657.74 1750 36
37 646.05 1800 37
38 706.03 1850 38
39 712.38 1900 39
40 718.78 1950 40
41 699.49 2000 41
42 635.36 2050 42
43 682.09 2100 43
44 722.70 2150 44
45 731.22 2200 45
46 763.95 2250 46
47 739.86 2300 47
48 744.88 2350 48
49 746.73 2400 49
50 821.77 2450 50
51 752.76 2500 51
52 733.80 2550 52
53 735.91 2600 53
54 783.64 2650 54
55 711.28 2700 55
56 764.41 2750 56
57 833.71 2800 57
58 827.09 2850 58
59 766.46 2900 59
60 748.42 2950 60
61 870.61 3000 61
62 854.52 3050 62
63 858.34 3100 63
64 787.99 3150 64
65 834.26 3200 65
66 827.86 3250 66
67 771.05 3300 67
68 806.20 3350 68
69 873.40 3400 69
70 792.56 3450 70
71 855.02 3500 71
72 794.63 3550 72
73 861.60 3600 73
74 859.60 3650 74
75 856.97 3700 75
76 905.18 3750 76
77 933.00 3800 77
78 838.89 3850 78
79 903.42 3900 79
80 889.50 3950 80
81 914.18 4000 81
82 863.96 4050 82
83 937.39 4000 83
84 948.76 4150 84
85 900.66 4200 85
86 947.49 4250 86
87 904.22 4300 87
88 861.64 4350 88
89 918.50 4400 89
90 906.68 4450 90
91 966.54 4500 91
92 997.92 4550 92
93 965.90 4600 93
94 969.30 4650 94
95 904.80 4700 95
96 957.80 4750 96
97 1026.98 4800 97
98 1048.42 4850 98
99 953.19 4900 99
100 1020.25 4950 100
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Mest t
480.3235 -0.2831 19.1481
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-79.240 -22.882 1.338 23.688 77.530
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 480.3235 18.8508 25.480 <2e-16 ***
Mest -0.2831 0.3562 -0.795 0.429
t 19.1481 17.7963 1.076 0.285
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 35.21 on 97 degrees of freedom
Multiple R-squared: 0.9456, Adjusted R-squared: 0.9444
F-statistic: 842.4 on 2 and 97 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.047389346 0.09477869 0.9526107
[2,] 0.018206736 0.03641347 0.9817933
[3,] 0.017672136 0.03534427 0.9823279
[4,] 0.011428816 0.02285763 0.9885712
[5,] 0.006646898 0.01329380 0.9933531
[6,] 0.008965525 0.01793105 0.9910345
[7,] 0.213044089 0.42608818 0.7869559
[8,] 0.144872101 0.28974420 0.8551279
[9,] 0.095294645 0.19058929 0.9047054
[10,] 0.213697071 0.42739414 0.7863029
[11,] 0.158169285 0.31633857 0.8418307
[12,] 0.109691720 0.21938344 0.8903083
[13,] 0.109656189 0.21931238 0.8903438
[14,] 0.094174855 0.18834971 0.9058251
[15,] 0.082573208 0.16514642 0.9174268
[16,] 0.075697135 0.15139427 0.9243029
[17,] 0.063225294 0.12645059 0.9367747
[18,] 0.045130804 0.09026161 0.9548692
[19,] 0.035783114 0.07156623 0.9642169
[20,] 0.024855262 0.04971052 0.9751447
[21,] 0.018148483 0.03629697 0.9818515
[22,] 0.017973933 0.03594787 0.9820261
[23,] 0.017231995 0.03446399 0.9827680
[24,] 0.015859214 0.03171843 0.9841408
[25,] 0.010566434 0.02113287 0.9894336
[26,] 0.017312390 0.03462478 0.9826876
[27,] 0.020690771 0.04138154 0.9793092
[28,] 0.164943603 0.32988721 0.8350564
[29,] 0.331361132 0.66272226 0.6686389
[30,] 0.280003500 0.56000700 0.7199965
[31,] 0.241997332 0.48399466 0.7580027
[32,] 0.236832804 0.47366561 0.7631672
[33,] 0.209210100 0.41842020 0.7907899
[34,] 0.183355526 0.36671105 0.8166445
[35,] 0.159825429 0.31965086 0.8401746
[36,] 0.125774949 0.25154990 0.8742251
[37,] 0.244125012 0.48825002 0.7558750
[38,] 0.227038916 0.45407783 0.7729611
[39,] 0.188251428 0.37650286 0.8117486
[40,] 0.155210261 0.31042052 0.8447897
[41,] 0.160675118 0.32135024 0.8393249
[42,] 0.128766124 0.25753225 0.8712339
[43,] 0.101427113 0.20285423 0.8985729
[44,] 0.077902060 0.15580412 0.9220979
[45,] 0.174983515 0.34996703 0.8250165
[46,] 0.140290729 0.28058146 0.8597093
[47,] 0.123887413 0.24777483 0.8761126
[48,] 0.111344957 0.22268991 0.8886550
[49,] 0.090683524 0.18136705 0.9093165
[50,] 0.145587207 0.29117441 0.8544128
[51,] 0.118494023 0.23698805 0.8815060
[52,] 0.150854032 0.30170806 0.8491460
[53,] 0.159513499 0.31902700 0.8404865
[54,] 0.141912590 0.28382518 0.8580874
[55,] 0.170350579 0.34070116 0.8296494
[56,] 0.288609177 0.57721835 0.7113908
[57,] 0.341753300 0.68350660 0.6582467
[58,] 0.412308966 0.82461793 0.5876910
[59,] 0.382115893 0.76423179 0.6178841
[60,] 0.349373647 0.69874729 0.6506264
[61,] 0.305711491 0.61142298 0.6942885
[62,] 0.372743749 0.74548750 0.6272563
[63,] 0.340630617 0.68126123 0.6593694
[64,] 0.357262501 0.71452500 0.6427375
[65,] 0.395172385 0.79034477 0.6048276
[66,] 0.341285907 0.68257181 0.6587141
[67,] 0.427440486 0.85488097 0.5725595
[68,] 0.363930111 0.72786022 0.6360699
[69,] 0.303085271 0.60617054 0.6969147
[70,] 0.251403679 0.50280736 0.7485963
[71,] 0.240554130 0.48110826 0.7594459
[72,] 0.352708803 0.70541761 0.6472912
[73,] 0.350836823 0.70167365 0.6491632
[74,] 0.312931213 0.62586243 0.6870688
[75,] 0.254164177 0.50832835 0.7458358
[76,] 0.235142639 0.47028528 0.7648574
[77,] 0.207140466 0.41428093 0.7928595
[78,] 0.155072628 0.31014526 0.8449274
[79,] 0.194314792 0.38862958 0.8056852
[80,] 0.144899857 0.28979971 0.8551001
[81,] 0.175710885 0.35142177 0.8242891
[82,] 0.128853240 0.25770648 0.8711468
[83,] 0.179739585 0.35947917 0.8202604
[84,] 0.127201748 0.25440350 0.8727983
[85,] 0.123003803 0.24600761 0.8769962
[86,] 0.078733347 0.15746669 0.9212667
[87,] 0.090903341 0.18180668 0.9090967
[88,] 0.058704378 0.11740876 0.9412956
[89,] 0.038523290 0.07704658 0.9614767
> postscript(file="/var/wessaorg/rcomp/tmp/14gkz1355319734.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/22zyv1355319734.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/3f4yc1355319734.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/4u23a1355319734.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/5baje1355319734.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 = 100
Frequency = 1
1 2 3 4 5
-4.132162e+01 -2.687689e+01 -4.552171e+00 -1.184745e+01 -5.332726e+00
6 7 8 9 10
-1.180800e+01 1.661672e+01 -9.178558e+00 3.176616e+01 6.790887e+00
11 12 13 14 15
5.483561e+01 -4.413967e+01 1.549506e+01 1.638978e+01 -4.207550e+01
16 17 18 19 20
-3.030776e+00 8.573946e+00 -2.674133e+01 3.022339e+01 3.672811e+01
21 22 23 24 25
-1.623716e+01 -1.497244e+01 -6.177718e+00 2.831701e+01 -3.448272e+00
26 27 28 29 30
2.620645e+01 -2.193883e+01 3.911590e+01 4.077062e+01 2.155534e+01
31 32 33 34 35
6.246006e+01 -1.803521e+01 -7.924049e+01 -6.867577e+01 -7.511045e+00
36 37 38 39 40
-1.656632e+01 -3.325160e+01 2.173312e+01 2.308785e+01 2.449257e+01
41 42 43 44 45
2.072909e-01 -6.891799e+01 -2.718326e+01 8.431459e+00 1.195618e+01
46 47 48 49 50
3.969090e+01 1.060563e+01 1.063035e+01 7.485072e+00 7.752980e+01
51 52 53 54 55
3.524518e+00 -2.043076e+01 -2.331604e+01 1.941869e+01 -5.793659e+01
56 57 58 59 60
-9.801869e+00 5.450285e+01 4.288758e+01 -2.273770e+01 -4.577298e+01
61 62 63 64 65
7.142174e+01 5.033647e+01 4.916119e+01 -2.618409e+01 1.509064e+01
66 67 68 69 70
3.695358e+00 -5.810992e+01 -2.795520e+01 3.424953e+01 -5.158575e+01
71 72 73 74 75
5.878972e+00 -5.950631e+01 2.468417e+00 -4.526860e+00 -1.215214e+01
76 77 78 79 80
3.106259e+01 5.388731e+01 -4.521797e+01 1.431675e+01 -4.598524e+00
81 82 83 84 85
1.508620e+01 -4.012908e+01 -6.360554e-13 3.468037e+01 -1.841491e+01
86 87 88 89 90
2.341981e+01 -2.484547e+01 -7.242074e+01 -2.055602e+01 -3.737130e+01
91 92 93 94 95
1.749343e+01 4.387815e+01 6.862871e+00 5.267594e+00 -6.422768e+01
96 97 98 99 100
-1.622296e+01 4.796176e+01 6.440648e+01 -3.581879e+01 2.624593e+01
> postscript(file="/var/wessaorg/rcomp/tmp/6dg7h1355319734.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 = 100
Frequency = 1
lag(myerror, k = 1) myerror
0 -4.132162e+01 NA
1 -2.687689e+01 -4.132162e+01
2 -4.552171e+00 -2.687689e+01
3 -1.184745e+01 -4.552171e+00
4 -5.332726e+00 -1.184745e+01
5 -1.180800e+01 -5.332726e+00
6 1.661672e+01 -1.180800e+01
7 -9.178558e+00 1.661672e+01
8 3.176616e+01 -9.178558e+00
9 6.790887e+00 3.176616e+01
10 5.483561e+01 6.790887e+00
11 -4.413967e+01 5.483561e+01
12 1.549506e+01 -4.413967e+01
13 1.638978e+01 1.549506e+01
14 -4.207550e+01 1.638978e+01
15 -3.030776e+00 -4.207550e+01
16 8.573946e+00 -3.030776e+00
17 -2.674133e+01 8.573946e+00
18 3.022339e+01 -2.674133e+01
19 3.672811e+01 3.022339e+01
20 -1.623716e+01 3.672811e+01
21 -1.497244e+01 -1.623716e+01
22 -6.177718e+00 -1.497244e+01
23 2.831701e+01 -6.177718e+00
24 -3.448272e+00 2.831701e+01
25 2.620645e+01 -3.448272e+00
26 -2.193883e+01 2.620645e+01
27 3.911590e+01 -2.193883e+01
28 4.077062e+01 3.911590e+01
29 2.155534e+01 4.077062e+01
30 6.246006e+01 2.155534e+01
31 -1.803521e+01 6.246006e+01
32 -7.924049e+01 -1.803521e+01
33 -6.867577e+01 -7.924049e+01
34 -7.511045e+00 -6.867577e+01
35 -1.656632e+01 -7.511045e+00
36 -3.325160e+01 -1.656632e+01
37 2.173312e+01 -3.325160e+01
38 2.308785e+01 2.173312e+01
39 2.449257e+01 2.308785e+01
40 2.072909e-01 2.449257e+01
41 -6.891799e+01 2.072909e-01
42 -2.718326e+01 -6.891799e+01
43 8.431459e+00 -2.718326e+01
44 1.195618e+01 8.431459e+00
45 3.969090e+01 1.195618e+01
46 1.060563e+01 3.969090e+01
47 1.063035e+01 1.060563e+01
48 7.485072e+00 1.063035e+01
49 7.752980e+01 7.485072e+00
50 3.524518e+00 7.752980e+01
51 -2.043076e+01 3.524518e+00
52 -2.331604e+01 -2.043076e+01
53 1.941869e+01 -2.331604e+01
54 -5.793659e+01 1.941869e+01
55 -9.801869e+00 -5.793659e+01
56 5.450285e+01 -9.801869e+00
57 4.288758e+01 5.450285e+01
58 -2.273770e+01 4.288758e+01
59 -4.577298e+01 -2.273770e+01
60 7.142174e+01 -4.577298e+01
61 5.033647e+01 7.142174e+01
62 4.916119e+01 5.033647e+01
63 -2.618409e+01 4.916119e+01
64 1.509064e+01 -2.618409e+01
65 3.695358e+00 1.509064e+01
66 -5.810992e+01 3.695358e+00
67 -2.795520e+01 -5.810992e+01
68 3.424953e+01 -2.795520e+01
69 -5.158575e+01 3.424953e+01
70 5.878972e+00 -5.158575e+01
71 -5.950631e+01 5.878972e+00
72 2.468417e+00 -5.950631e+01
73 -4.526860e+00 2.468417e+00
74 -1.215214e+01 -4.526860e+00
75 3.106259e+01 -1.215214e+01
76 5.388731e+01 3.106259e+01
77 -4.521797e+01 5.388731e+01
78 1.431675e+01 -4.521797e+01
79 -4.598524e+00 1.431675e+01
80 1.508620e+01 -4.598524e+00
81 -4.012908e+01 1.508620e+01
82 -6.360554e-13 -4.012908e+01
83 3.468037e+01 -6.360554e-13
84 -1.841491e+01 3.468037e+01
85 2.341981e+01 -1.841491e+01
86 -2.484547e+01 2.341981e+01
87 -7.242074e+01 -2.484547e+01
88 -2.055602e+01 -7.242074e+01
89 -3.737130e+01 -2.055602e+01
90 1.749343e+01 -3.737130e+01
91 4.387815e+01 1.749343e+01
92 6.862871e+00 4.387815e+01
93 5.267594e+00 6.862871e+00
94 -6.422768e+01 5.267594e+00
95 -1.622296e+01 -6.422768e+01
96 4.796176e+01 -1.622296e+01
97 6.440648e+01 4.796176e+01
98 -3.581879e+01 6.440648e+01
99 2.624593e+01 -3.581879e+01
100 NA 2.624593e+01
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -2.687689e+01 -4.132162e+01
[2,] -4.552171e+00 -2.687689e+01
[3,] -1.184745e+01 -4.552171e+00
[4,] -5.332726e+00 -1.184745e+01
[5,] -1.180800e+01 -5.332726e+00
[6,] 1.661672e+01 -1.180800e+01
[7,] -9.178558e+00 1.661672e+01
[8,] 3.176616e+01 -9.178558e+00
[9,] 6.790887e+00 3.176616e+01
[10,] 5.483561e+01 6.790887e+00
[11,] -4.413967e+01 5.483561e+01
[12,] 1.549506e+01 -4.413967e+01
[13,] 1.638978e+01 1.549506e+01
[14,] -4.207550e+01 1.638978e+01
[15,] -3.030776e+00 -4.207550e+01
[16,] 8.573946e+00 -3.030776e+00
[17,] -2.674133e+01 8.573946e+00
[18,] 3.022339e+01 -2.674133e+01
[19,] 3.672811e+01 3.022339e+01
[20,] -1.623716e+01 3.672811e+01
[21,] -1.497244e+01 -1.623716e+01
[22,] -6.177718e+00 -1.497244e+01
[23,] 2.831701e+01 -6.177718e+00
[24,] -3.448272e+00 2.831701e+01
[25,] 2.620645e+01 -3.448272e+00
[26,] -2.193883e+01 2.620645e+01
[27,] 3.911590e+01 -2.193883e+01
[28,] 4.077062e+01 3.911590e+01
[29,] 2.155534e+01 4.077062e+01
[30,] 6.246006e+01 2.155534e+01
[31,] -1.803521e+01 6.246006e+01
[32,] -7.924049e+01 -1.803521e+01
[33,] -6.867577e+01 -7.924049e+01
[34,] -7.511045e+00 -6.867577e+01
[35,] -1.656632e+01 -7.511045e+00
[36,] -3.325160e+01 -1.656632e+01
[37,] 2.173312e+01 -3.325160e+01
[38,] 2.308785e+01 2.173312e+01
[39,] 2.449257e+01 2.308785e+01
[40,] 2.072909e-01 2.449257e+01
[41,] -6.891799e+01 2.072909e-01
[42,] -2.718326e+01 -6.891799e+01
[43,] 8.431459e+00 -2.718326e+01
[44,] 1.195618e+01 8.431459e+00
[45,] 3.969090e+01 1.195618e+01
[46,] 1.060563e+01 3.969090e+01
[47,] 1.063035e+01 1.060563e+01
[48,] 7.485072e+00 1.063035e+01
[49,] 7.752980e+01 7.485072e+00
[50,] 3.524518e+00 7.752980e+01
[51,] -2.043076e+01 3.524518e+00
[52,] -2.331604e+01 -2.043076e+01
[53,] 1.941869e+01 -2.331604e+01
[54,] -5.793659e+01 1.941869e+01
[55,] -9.801869e+00 -5.793659e+01
[56,] 5.450285e+01 -9.801869e+00
[57,] 4.288758e+01 5.450285e+01
[58,] -2.273770e+01 4.288758e+01
[59,] -4.577298e+01 -2.273770e+01
[60,] 7.142174e+01 -4.577298e+01
[61,] 5.033647e+01 7.142174e+01
[62,] 4.916119e+01 5.033647e+01
[63,] -2.618409e+01 4.916119e+01
[64,] 1.509064e+01 -2.618409e+01
[65,] 3.695358e+00 1.509064e+01
[66,] -5.810992e+01 3.695358e+00
[67,] -2.795520e+01 -5.810992e+01
[68,] 3.424953e+01 -2.795520e+01
[69,] -5.158575e+01 3.424953e+01
[70,] 5.878972e+00 -5.158575e+01
[71,] -5.950631e+01 5.878972e+00
[72,] 2.468417e+00 -5.950631e+01
[73,] -4.526860e+00 2.468417e+00
[74,] -1.215214e+01 -4.526860e+00
[75,] 3.106259e+01 -1.215214e+01
[76,] 5.388731e+01 3.106259e+01
[77,] -4.521797e+01 5.388731e+01
[78,] 1.431675e+01 -4.521797e+01
[79,] -4.598524e+00 1.431675e+01
[80,] 1.508620e+01 -4.598524e+00
[81,] -4.012908e+01 1.508620e+01
[82,] -6.360554e-13 -4.012908e+01
[83,] 3.468037e+01 -6.360554e-13
[84,] -1.841491e+01 3.468037e+01
[85,] 2.341981e+01 -1.841491e+01
[86,] -2.484547e+01 2.341981e+01
[87,] -7.242074e+01 -2.484547e+01
[88,] -2.055602e+01 -7.242074e+01
[89,] -3.737130e+01 -2.055602e+01
[90,] 1.749343e+01 -3.737130e+01
[91,] 4.387815e+01 1.749343e+01
[92,] 6.862871e+00 4.387815e+01
[93,] 5.267594e+00 6.862871e+00
[94,] -6.422768e+01 5.267594e+00
[95,] -1.622296e+01 -6.422768e+01
[96,] 4.796176e+01 -1.622296e+01
[97,] 6.440648e+01 4.796176e+01
[98,] -3.581879e+01 6.440648e+01
[99,] 2.624593e+01 -3.581879e+01
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -2.687689e+01 -4.132162e+01
2 -4.552171e+00 -2.687689e+01
3 -1.184745e+01 -4.552171e+00
4 -5.332726e+00 -1.184745e+01
5 -1.180800e+01 -5.332726e+00
6 1.661672e+01 -1.180800e+01
7 -9.178558e+00 1.661672e+01
8 3.176616e+01 -9.178558e+00
9 6.790887e+00 3.176616e+01
10 5.483561e+01 6.790887e+00
11 -4.413967e+01 5.483561e+01
12 1.549506e+01 -4.413967e+01
13 1.638978e+01 1.549506e+01
14 -4.207550e+01 1.638978e+01
15 -3.030776e+00 -4.207550e+01
16 8.573946e+00 -3.030776e+00
17 -2.674133e+01 8.573946e+00
18 3.022339e+01 -2.674133e+01
19 3.672811e+01 3.022339e+01
20 -1.623716e+01 3.672811e+01
21 -1.497244e+01 -1.623716e+01
22 -6.177718e+00 -1.497244e+01
23 2.831701e+01 -6.177718e+00
24 -3.448272e+00 2.831701e+01
25 2.620645e+01 -3.448272e+00
26 -2.193883e+01 2.620645e+01
27 3.911590e+01 -2.193883e+01
28 4.077062e+01 3.911590e+01
29 2.155534e+01 4.077062e+01
30 6.246006e+01 2.155534e+01
31 -1.803521e+01 6.246006e+01
32 -7.924049e+01 -1.803521e+01
33 -6.867577e+01 -7.924049e+01
34 -7.511045e+00 -6.867577e+01
35 -1.656632e+01 -7.511045e+00
36 -3.325160e+01 -1.656632e+01
37 2.173312e+01 -3.325160e+01
38 2.308785e+01 2.173312e+01
39 2.449257e+01 2.308785e+01
40 2.072909e-01 2.449257e+01
41 -6.891799e+01 2.072909e-01
42 -2.718326e+01 -6.891799e+01
43 8.431459e+00 -2.718326e+01
44 1.195618e+01 8.431459e+00
45 3.969090e+01 1.195618e+01
46 1.060563e+01 3.969090e+01
47 1.063035e+01 1.060563e+01
48 7.485072e+00 1.063035e+01
49 7.752980e+01 7.485072e+00
50 3.524518e+00 7.752980e+01
51 -2.043076e+01 3.524518e+00
52 -2.331604e+01 -2.043076e+01
53 1.941869e+01 -2.331604e+01
54 -5.793659e+01 1.941869e+01
55 -9.801869e+00 -5.793659e+01
56 5.450285e+01 -9.801869e+00
57 4.288758e+01 5.450285e+01
58 -2.273770e+01 4.288758e+01
59 -4.577298e+01 -2.273770e+01
60 7.142174e+01 -4.577298e+01
61 5.033647e+01 7.142174e+01
62 4.916119e+01 5.033647e+01
63 -2.618409e+01 4.916119e+01
64 1.509064e+01 -2.618409e+01
65 3.695358e+00 1.509064e+01
66 -5.810992e+01 3.695358e+00
67 -2.795520e+01 -5.810992e+01
68 3.424953e+01 -2.795520e+01
69 -5.158575e+01 3.424953e+01
70 5.878972e+00 -5.158575e+01
71 -5.950631e+01 5.878972e+00
72 2.468417e+00 -5.950631e+01
73 -4.526860e+00 2.468417e+00
74 -1.215214e+01 -4.526860e+00
75 3.106259e+01 -1.215214e+01
76 5.388731e+01 3.106259e+01
77 -4.521797e+01 5.388731e+01
78 1.431675e+01 -4.521797e+01
79 -4.598524e+00 1.431675e+01
80 1.508620e+01 -4.598524e+00
81 -4.012908e+01 1.508620e+01
82 -6.360554e-13 -4.012908e+01
83 3.468037e+01 -6.360554e-13
84 -1.841491e+01 3.468037e+01
85 2.341981e+01 -1.841491e+01
86 -2.484547e+01 2.341981e+01
87 -7.242074e+01 -2.484547e+01
88 -2.055602e+01 -7.242074e+01
89 -3.737130e+01 -2.055602e+01
90 1.749343e+01 -3.737130e+01
91 4.387815e+01 1.749343e+01
92 6.862871e+00 4.387815e+01
93 5.267594e+00 6.862871e+00
94 -6.422768e+01 5.267594e+00
95 -1.622296e+01 -6.422768e+01
96 4.796176e+01 -1.622296e+01
97 6.440648e+01 4.796176e+01
98 -3.581879e+01 6.440648e+01
99 2.624593e+01 -3.581879e+01
> 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/7egrf1355319734.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/8tuz91355319734.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/9juyd1355319734.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')
Warning messages:
1: Not plotting observations with leverage one:
83
2: Not plotting observations with leverage one:
83
> par(opar)
> dev.off()
null device
1
> if (n > n25) {
+ postscript(file="/var/wessaorg/rcomp/tmp/1041pn1355319734.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/112ce31355319734.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/12qhj91355319734.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/139sgm1355319734.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/14at1v1355319734.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/15q9g01355319734.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/16jtyr1355319734.tab")
+ }
>
> try(system("convert tmp/14gkz1355319734.ps tmp/14gkz1355319734.png",intern=TRUE))
character(0)
> try(system("convert tmp/22zyv1355319734.ps tmp/22zyv1355319734.png",intern=TRUE))
character(0)
> try(system("convert tmp/3f4yc1355319734.ps tmp/3f4yc1355319734.png",intern=TRUE))
character(0)
> try(system("convert tmp/4u23a1355319734.ps tmp/4u23a1355319734.png",intern=TRUE))
character(0)
> try(system("convert tmp/5baje1355319734.ps tmp/5baje1355319734.png",intern=TRUE))
character(0)
> try(system("convert tmp/6dg7h1355319734.ps tmp/6dg7h1355319734.png",intern=TRUE))
character(0)
> try(system("convert tmp/7egrf1355319734.ps tmp/7egrf1355319734.png",intern=TRUE))
character(0)
> try(system("convert tmp/8tuz91355319734.ps tmp/8tuz91355319734.png",intern=TRUE))
character(0)
> try(system("convert tmp/9juyd1355319734.ps tmp/9juyd1355319734.png",intern=TRUE))
character(0)
> try(system("convert tmp/1041pn1355319734.ps tmp/1041pn1355319734.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
6.917 1.117 8.026