R version 2.8.0 (2008-10-20)
Copyright (C) 2008 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> x <- array(list(11554.5
+ ,180144
+ ,13182.1
+ ,173666
+ ,14800.1
+ ,165688
+ ,12150.7
+ ,161570
+ ,14478.2
+ ,156145
+ ,13253.9
+ ,153730
+ ,12036.8
+ ,182698
+ ,12653.2
+ ,200765
+ ,14035.4
+ ,176512
+ ,14571.4
+ ,166618
+ ,15400.9
+ ,158644
+ ,14283.2
+ ,159585
+ ,14485.3
+ ,163095
+ ,14196.3
+ ,159044
+ ,15559.1
+ ,155511
+ ,13767.4
+ ,153745
+ ,14634
+ ,150569
+ ,14381.1
+ ,150605
+ ,12509.9
+ ,179612
+ ,12122.3
+ ,194690
+ ,13122.3
+ ,189917
+ ,13908.7
+ ,184128
+ ,13456.5
+ ,175335
+ ,12441.6
+ ,179566
+ ,12953
+ ,181140
+ ,13057.2
+ ,177876
+ ,14350.1
+ ,175041
+ ,13830.2
+ ,169292
+ ,13755.5
+ ,166070
+ ,13574.4
+ ,166972
+ ,12802.6
+ ,206348
+ ,11737.3
+ ,215706
+ ,13850.2
+ ,202108
+ ,15081.8
+ ,195411
+ ,13653.3
+ ,193111
+ ,14019.1
+ ,195198
+ ,13962
+ ,198770
+ ,13768.7
+ ,194163
+ ,14747.1
+ ,190420
+ ,13858.1
+ ,189733
+ ,13188
+ ,186029
+ ,13693.1
+ ,191531
+ ,12970
+ ,232571
+ ,11392.8
+ ,243477
+ ,13985.2
+ ,227247
+ ,14994.7
+ ,217859
+ ,13584.7
+ ,208679
+ ,14257.8
+ ,213188
+ ,13553.4
+ ,216234
+ ,14007.3
+ ,213586
+ ,16535.8
+ ,209465
+ ,14721.4
+ ,204045
+ ,13664.6
+ ,200237
+ ,16805.9
+ ,203666
+ ,13829.4
+ ,241476
+ ,13735.6
+ ,260307
+ ,15870.5
+ ,243324
+ ,15962.4
+ ,244460
+ ,15744.1
+ ,233575
+ ,16083.7
+ ,237217
+ ,14863.9
+ ,235243
+ ,15533.1
+ ,230354
+ ,17473.1
+ ,227184
+ ,15925.5
+ ,221678
+ ,15573.7
+ ,217142
+ ,17495
+ ,219452
+ ,14155.8
+ ,256446
+ ,14913.9
+ ,265845
+ ,17250.4
+ ,248624
+ ,15879.8
+ ,241114
+ ,17647.8
+ ,229245
+ ,17749.9
+ ,231805
+ ,17111.8
+ ,219277
+ ,16934.8
+ ,219313
+ ,20280
+ ,212610
+ ,16238.2
+ ,214771
+ ,17896.1
+ ,211142
+ ,18089.3
+ ,211457
+ ,15660
+ ,240048
+ ,16162.4
+ ,240636
+ ,17850.1
+ ,230580
+ ,18520.4
+ ,208795
+ ,18524.7
+ ,197922
+ ,16843.7
+ ,194596)
+ ,dim=c(2
+ ,84)
+ ,dimnames=list(c('invoer'
+ ,'werkloosheid')
+ ,1:84))
> y <- array(NA,dim=c(2,84),dimnames=list(c('invoer','werkloosheid'),1:84))
> 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
Attaching package: 'zoo'
The following object(s) are masked from package:base :
as.Date.numeric
> n25 <- 25 #minimum number of obs. for Goldfeld-Quandt test
> par1 <- as.numeric(par1)
> x <- t(y)
> k <- length(x[1,])
> n <- length(x[,1])
> x1 <- cbind(x[,par1], x[,1:k!=par1])
> mycolnames <- c(colnames(x)[par1], colnames(x)[1:k!=par1])
> colnames(x1) <- mycolnames #colnames(x)[par1]
> x <- x1
> if (par3 == 'First Differences'){
+ x2 <- array(0, dim=c(n-1,k), dimnames=list(1:(n-1), paste('(1-B)',colnames(x),sep='')))
+ for (i in 1:n-1) {
+ for (j in 1:k) {
+ x2[i,j] <- x[i+1,j] - x[i,j]
+ }
+ }
+ x <- x2
+ }
> if (par2 == 'Include Monthly Dummies'){
+ x2 <- array(0, dim=c(n,11), dimnames=list(1:n, paste('M', seq(1:11), sep ='')))
+ for (i in 1:11){
+ x2[seq(i,n,12),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> if (par2 == 'Include Quarterly Dummies'){
+ x2 <- array(0, dim=c(n,3), dimnames=list(1:n, paste('Q', seq(1:3), sep ='')))
+ for (i in 1:3){
+ x2[seq(i,n,4),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> k <- length(x[1,])
> if (par3 == 'Linear Trend'){
+ x <- cbind(x, c(1:n))
+ colnames(x)[k+1] <- 't'
+ }
> x
invoer werkloosheid
1 11554.5 180144
2 13182.1 173666
3 14800.1 165688
4 12150.7 161570
5 14478.2 156145
6 13253.9 153730
7 12036.8 182698
8 12653.2 200765
9 14035.4 176512
10 14571.4 166618
11 15400.9 158644
12 14283.2 159585
13 14485.3 163095
14 14196.3 159044
15 15559.1 155511
16 13767.4 153745
17 14634.0 150569
18 14381.1 150605
19 12509.9 179612
20 12122.3 194690
21 13122.3 189917
22 13908.7 184128
23 13456.5 175335
24 12441.6 179566
25 12953.0 181140
26 13057.2 177876
27 14350.1 175041
28 13830.2 169292
29 13755.5 166070
30 13574.4 166972
31 12802.6 206348
32 11737.3 215706
33 13850.2 202108
34 15081.8 195411
35 13653.3 193111
36 14019.1 195198
37 13962.0 198770
38 13768.7 194163
39 14747.1 190420
40 13858.1 189733
41 13188.0 186029
42 13693.1 191531
43 12970.0 232571
44 11392.8 243477
45 13985.2 227247
46 14994.7 217859
47 13584.7 208679
48 14257.8 213188
49 13553.4 216234
50 14007.3 213586
51 16535.8 209465
52 14721.4 204045
53 13664.6 200237
54 16805.9 203666
55 13829.4 241476
56 13735.6 260307
57 15870.5 243324
58 15962.4 244460
59 15744.1 233575
60 16083.7 237217
61 14863.9 235243
62 15533.1 230354
63 17473.1 227184
64 15925.5 221678
65 15573.7 217142
66 17495.0 219452
67 14155.8 256446
68 14913.9 265845
69 17250.4 248624
70 15879.8 241114
71 17647.8 229245
72 17749.9 231805
73 17111.8 219277
74 16934.8 219313
75 20280.0 212610
76 16238.2 214771
77 17896.1 211142
78 18089.3 211457
79 15660.0 240048
80 16162.4 240636
81 17850.1 230580
82 18520.4 208795
83 18524.7 197922
84 16843.7 194596
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) werkloosheid
1.034e+04 2.187e-02
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-4272.46 -1185.50 -96.02 768.88 5289.68
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.034e+04 1.307e+03 7.914 1.03e-11 ***
werkloosheid 2.187e-02 6.411e-03 3.410 0.00101 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1720 on 82 degrees of freedom
Multiple R-squared: 0.1242, Adjusted R-squared: 0.1135
F-statistic: 11.63 on 1 and 82 DF, p-value: 0.001009
> 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,] 3.689488e-01 7.378975e-01 0.63105124
[2,] 2.616203e-01 5.232407e-01 0.73837967
[3,] 1.544979e-01 3.089958e-01 0.84550209
[4,] 1.254241e-01 2.508481e-01 0.87457593
[5,] 1.010291e-01 2.020582e-01 0.89897091
[6,] 8.495495e-02 1.699099e-01 0.91504505
[7,] 9.113162e-02 1.822632e-01 0.90886838
[8,] 5.439982e-02 1.087996e-01 0.94560018
[9,] 3.408932e-02 6.817864e-02 0.96591068
[10,] 1.854661e-02 3.709321e-02 0.98145339
[11,] 1.662473e-02 3.324946e-02 0.98337527
[12,] 1.062002e-02 2.124004e-02 0.98937998
[13,] 5.664868e-03 1.132974e-02 0.99433513
[14,] 2.965797e-03 5.931594e-03 0.99703420
[15,] 1.788430e-03 3.576859e-03 0.99821157
[16,] 1.051652e-03 2.103304e-03 0.99894835
[17,] 6.348848e-04 1.269770e-03 0.99936512
[18,] 4.777720e-04 9.555441e-04 0.99952223
[19,] 2.286938e-04 4.573876e-04 0.99977131
[20,] 1.585306e-04 3.170611e-04 0.99984147
[21,] 8.101934e-05 1.620387e-04 0.99991898
[22,] 4.075368e-05 8.150735e-05 0.99995925
[23,] 2.998280e-05 5.996559e-05 0.99997002
[24,] 1.356209e-05 2.712418e-05 0.99998644
[25,] 6.013796e-06 1.202759e-05 0.99999399
[26,] 2.793525e-06 5.587049e-06 0.99999721
[27,] 2.543239e-06 5.086478e-06 0.99999746
[28,] 2.566332e-06 5.132664e-06 0.99999743
[29,] 4.638783e-06 9.277565e-06 0.99999536
[30,] 2.541272e-05 5.082544e-05 0.99997459
[31,] 1.812974e-05 3.625948e-05 0.99998187
[32,] 1.600006e-05 3.200012e-05 0.99998400
[33,] 1.431419e-05 2.862838e-05 0.99998569
[34,] 1.095664e-05 2.191329e-05 0.99998904
[35,] 1.374229e-05 2.748459e-05 0.99998626
[36,] 1.107984e-05 2.215969e-05 0.99998892
[37,] 1.428404e-05 2.856808e-05 0.99998572
[38,] 1.802698e-05 3.605396e-05 0.99998197
[39,] 2.141540e-05 4.283079e-05 0.99997858
[40,] 8.935478e-05 1.787096e-04 0.99991065
[41,] 1.614260e-04 3.228521e-04 0.99983857
[42,] 4.065507e-04 8.131014e-04 0.99959345
[43,] 7.253642e-04 1.450728e-03 0.99927464
[44,] 1.178088e-03 2.356176e-03 0.99882191
[45,] 2.568995e-03 5.137990e-03 0.99743100
[46,] 5.619899e-03 1.123980e-02 0.99438010
[47,] 2.302699e-02 4.605398e-02 0.97697301
[48,] 4.565333e-02 9.130666e-02 0.95434667
[49,] 2.649037e-01 5.298074e-01 0.73509628
[50,] 4.482561e-01 8.965122e-01 0.55174389
[51,] 5.365211e-01 9.269578e-01 0.46347889
[52,] 5.324410e-01 9.351179e-01 0.46755896
[53,] 5.583390e-01 8.833220e-01 0.44166099
[54,] 5.667109e-01 8.665781e-01 0.43328905
[55,] 5.660851e-01 8.678298e-01 0.43391490
[56,] 5.552028e-01 8.895944e-01 0.44479720
[57,] 5.886142e-01 8.227717e-01 0.41138583
[58,] 6.012163e-01 7.975675e-01 0.39878373
[59,] 6.597906e-01 6.804189e-01 0.34020945
[60,] 6.784283e-01 6.431434e-01 0.32157171
[61,] 7.702321e-01 4.595357e-01 0.22976785
[62,] 7.764529e-01 4.470942e-01 0.22354710
[63,] 8.109074e-01 3.781851e-01 0.18909257
[64,] 7.719343e-01 4.561315e-01 0.22806574
[65,] 7.722632e-01 4.554736e-01 0.22773681
[66,] 7.312526e-01 5.374948e-01 0.26874738
[67,] 7.135315e-01 5.729371e-01 0.28646854
[68,] 7.045291e-01 5.909419e-01 0.29547095
[69,] 6.447712e-01 7.104576e-01 0.35522880
[70,] 5.782579e-01 8.434843e-01 0.42174213
[71,] 9.151032e-01 1.697935e-01 0.08489676
[72,] 9.163987e-01 1.672026e-01 0.08360129
[73,] 8.605635e-01 2.788730e-01 0.13943648
[74,] 7.862207e-01 4.275586e-01 0.21377930
[75,] 7.227314e-01 5.545373e-01 0.27726863
> postscript(file="/var/www/html/freestat/rcomp/tmp/1sauh1229261596.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> plot(x[,1], type='l', main='Actuals and Interpolation', ylab='value of Actuals and Interpolation (dots)', xlab='time or index')
> points(x[,1]-mysum$resid)
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/2o1hb1229261596.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> plot(mysum$resid, type='b', pch=19, main='Residuals', ylab='value of Residuals', xlab='time or index')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/3edx11229261596.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> hist(mysum$resid, main='Residual Histogram', xlab='values of Residuals')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/4h1d51229261596.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> densityplot(~mysum$resid,col='black',main='Residual Density Plot', xlab='values of Residuals')
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/5nno81229261596.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> qqnorm(mysum$resid, main='Residual Normal Q-Q Plot')
> qqline(mysum$resid)
> grid()
> dev.off()
null device
1
> (myerror <- as.ts(mysum$resid))
Time Series:
Start = 1
End = 84
Frequency = 1
1 2 3 4 5 6
-2725.92011 -956.67195 835.77524 -1723.58044 722.54278 -448.95075
7 8 9 10 11 12
-2299.46595 -2078.11954 -165.60268 586.73983 1590.59957 452.32363
13 14 15 16 17 18
577.67386 377.25316 1817.30585 64.22126 1000.26777 746.58060
19 20 21 22 23 24
-1758.88738 -2476.18340 -1371.81683 -458.83438 -718.76636 -1826.18154
25 26 27 28 29 30
-1349.19867 -1173.62795 181.26225 -212.92993 -217.17758 -418.00074
31 32 33 34 35 36
-2050.79759 -3320.71992 -910.48562 467.55119 -910.65694 -590.49134
37 38 39 40 41 42
-725.69681 -818.26000 241.98455 -631.99348 -1221.10170 -836.30861
43 44 45 46 47 48
-2456.79054 -4272.46147 -1325.17578 -110.39748 -1319.66731 -745.16125
49 50 51 52 53 54
-1516.16518 -1004.36393 1614.24599 -81.64012 -1055.17428 2011.14710
55 56 57 58 59 60
-1792.10754 -2297.66678 208.58403 275.64422 295.35597 555.31988
61 62 63 64 65 66
-621.31657 154.78646 2164.10178 736.89615 484.28046 2355.06993
67 68 69 70 71 72
-1793.04204 -1240.46086 1472.59406 266.20796 2293.73589 2339.85884
73 74 75 76 77 78
1975.69649 1797.90931 5289.67732 1200.62482 2937.87665 3124.18885
79 80 81 82 83 84
69.71715 559.25992 2466.84473 3613.49623 3855.54559 2247.27201
> postscript(file="/var/www/html/freestat/rcomp/tmp/65cgk1229261596.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> dum <- cbind(lag(myerror,k=1),myerror)
> dum
Time Series:
Start = 0
End = 84
Frequency = 1
lag(myerror, k = 1) myerror
0 -2725.92011 NA
1 -956.67195 -2725.92011
2 835.77524 -956.67195
3 -1723.58044 835.77524
4 722.54278 -1723.58044
5 -448.95075 722.54278
6 -2299.46595 -448.95075
7 -2078.11954 -2299.46595
8 -165.60268 -2078.11954
9 586.73983 -165.60268
10 1590.59957 586.73983
11 452.32363 1590.59957
12 577.67386 452.32363
13 377.25316 577.67386
14 1817.30585 377.25316
15 64.22126 1817.30585
16 1000.26777 64.22126
17 746.58060 1000.26777
18 -1758.88738 746.58060
19 -2476.18340 -1758.88738
20 -1371.81683 -2476.18340
21 -458.83438 -1371.81683
22 -718.76636 -458.83438
23 -1826.18154 -718.76636
24 -1349.19867 -1826.18154
25 -1173.62795 -1349.19867
26 181.26225 -1173.62795
27 -212.92993 181.26225
28 -217.17758 -212.92993
29 -418.00074 -217.17758
30 -2050.79759 -418.00074
31 -3320.71992 -2050.79759
32 -910.48562 -3320.71992
33 467.55119 -910.48562
34 -910.65694 467.55119
35 -590.49134 -910.65694
36 -725.69681 -590.49134
37 -818.26000 -725.69681
38 241.98455 -818.26000
39 -631.99348 241.98455
40 -1221.10170 -631.99348
41 -836.30861 -1221.10170
42 -2456.79054 -836.30861
43 -4272.46147 -2456.79054
44 -1325.17578 -4272.46147
45 -110.39748 -1325.17578
46 -1319.66731 -110.39748
47 -745.16125 -1319.66731
48 -1516.16518 -745.16125
49 -1004.36393 -1516.16518
50 1614.24599 -1004.36393
51 -81.64012 1614.24599
52 -1055.17428 -81.64012
53 2011.14710 -1055.17428
54 -1792.10754 2011.14710
55 -2297.66678 -1792.10754
56 208.58403 -2297.66678
57 275.64422 208.58403
58 295.35597 275.64422
59 555.31988 295.35597
60 -621.31657 555.31988
61 154.78646 -621.31657
62 2164.10178 154.78646
63 736.89615 2164.10178
64 484.28046 736.89615
65 2355.06993 484.28046
66 -1793.04204 2355.06993
67 -1240.46086 -1793.04204
68 1472.59406 -1240.46086
69 266.20796 1472.59406
70 2293.73589 266.20796
71 2339.85884 2293.73589
72 1975.69649 2339.85884
73 1797.90931 1975.69649
74 5289.67732 1797.90931
75 1200.62482 5289.67732
76 2937.87665 1200.62482
77 3124.18885 2937.87665
78 69.71715 3124.18885
79 559.25992 69.71715
80 2466.84473 559.25992
81 3613.49623 2466.84473
82 3855.54559 3613.49623
83 2247.27201 3855.54559
84 NA 2247.27201
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -956.67195 -2725.92011
[2,] 835.77524 -956.67195
[3,] -1723.58044 835.77524
[4,] 722.54278 -1723.58044
[5,] -448.95075 722.54278
[6,] -2299.46595 -448.95075
[7,] -2078.11954 -2299.46595
[8,] -165.60268 -2078.11954
[9,] 586.73983 -165.60268
[10,] 1590.59957 586.73983
[11,] 452.32363 1590.59957
[12,] 577.67386 452.32363
[13,] 377.25316 577.67386
[14,] 1817.30585 377.25316
[15,] 64.22126 1817.30585
[16,] 1000.26777 64.22126
[17,] 746.58060 1000.26777
[18,] -1758.88738 746.58060
[19,] -2476.18340 -1758.88738
[20,] -1371.81683 -2476.18340
[21,] -458.83438 -1371.81683
[22,] -718.76636 -458.83438
[23,] -1826.18154 -718.76636
[24,] -1349.19867 -1826.18154
[25,] -1173.62795 -1349.19867
[26,] 181.26225 -1173.62795
[27,] -212.92993 181.26225
[28,] -217.17758 -212.92993
[29,] -418.00074 -217.17758
[30,] -2050.79759 -418.00074
[31,] -3320.71992 -2050.79759
[32,] -910.48562 -3320.71992
[33,] 467.55119 -910.48562
[34,] -910.65694 467.55119
[35,] -590.49134 -910.65694
[36,] -725.69681 -590.49134
[37,] -818.26000 -725.69681
[38,] 241.98455 -818.26000
[39,] -631.99348 241.98455
[40,] -1221.10170 -631.99348
[41,] -836.30861 -1221.10170
[42,] -2456.79054 -836.30861
[43,] -4272.46147 -2456.79054
[44,] -1325.17578 -4272.46147
[45,] -110.39748 -1325.17578
[46,] -1319.66731 -110.39748
[47,] -745.16125 -1319.66731
[48,] -1516.16518 -745.16125
[49,] -1004.36393 -1516.16518
[50,] 1614.24599 -1004.36393
[51,] -81.64012 1614.24599
[52,] -1055.17428 -81.64012
[53,] 2011.14710 -1055.17428
[54,] -1792.10754 2011.14710
[55,] -2297.66678 -1792.10754
[56,] 208.58403 -2297.66678
[57,] 275.64422 208.58403
[58,] 295.35597 275.64422
[59,] 555.31988 295.35597
[60,] -621.31657 555.31988
[61,] 154.78646 -621.31657
[62,] 2164.10178 154.78646
[63,] 736.89615 2164.10178
[64,] 484.28046 736.89615
[65,] 2355.06993 484.28046
[66,] -1793.04204 2355.06993
[67,] -1240.46086 -1793.04204
[68,] 1472.59406 -1240.46086
[69,] 266.20796 1472.59406
[70,] 2293.73589 266.20796
[71,] 2339.85884 2293.73589
[72,] 1975.69649 2339.85884
[73,] 1797.90931 1975.69649
[74,] 5289.67732 1797.90931
[75,] 1200.62482 5289.67732
[76,] 2937.87665 1200.62482
[77,] 3124.18885 2937.87665
[78,] 69.71715 3124.18885
[79,] 559.25992 69.71715
[80,] 2466.84473 559.25992
[81,] 3613.49623 2466.84473
[82,] 3855.54559 3613.49623
[83,] 2247.27201 3855.54559
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -956.67195 -2725.92011
2 835.77524 -956.67195
3 -1723.58044 835.77524
4 722.54278 -1723.58044
5 -448.95075 722.54278
6 -2299.46595 -448.95075
7 -2078.11954 -2299.46595
8 -165.60268 -2078.11954
9 586.73983 -165.60268
10 1590.59957 586.73983
11 452.32363 1590.59957
12 577.67386 452.32363
13 377.25316 577.67386
14 1817.30585 377.25316
15 64.22126 1817.30585
16 1000.26777 64.22126
17 746.58060 1000.26777
18 -1758.88738 746.58060
19 -2476.18340 -1758.88738
20 -1371.81683 -2476.18340
21 -458.83438 -1371.81683
22 -718.76636 -458.83438
23 -1826.18154 -718.76636
24 -1349.19867 -1826.18154
25 -1173.62795 -1349.19867
26 181.26225 -1173.62795
27 -212.92993 181.26225
28 -217.17758 -212.92993
29 -418.00074 -217.17758
30 -2050.79759 -418.00074
31 -3320.71992 -2050.79759
32 -910.48562 -3320.71992
33 467.55119 -910.48562
34 -910.65694 467.55119
35 -590.49134 -910.65694
36 -725.69681 -590.49134
37 -818.26000 -725.69681
38 241.98455 -818.26000
39 -631.99348 241.98455
40 -1221.10170 -631.99348
41 -836.30861 -1221.10170
42 -2456.79054 -836.30861
43 -4272.46147 -2456.79054
44 -1325.17578 -4272.46147
45 -110.39748 -1325.17578
46 -1319.66731 -110.39748
47 -745.16125 -1319.66731
48 -1516.16518 -745.16125
49 -1004.36393 -1516.16518
50 1614.24599 -1004.36393
51 -81.64012 1614.24599
52 -1055.17428 -81.64012
53 2011.14710 -1055.17428
54 -1792.10754 2011.14710
55 -2297.66678 -1792.10754
56 208.58403 -2297.66678
57 275.64422 208.58403
58 295.35597 275.64422
59 555.31988 295.35597
60 -621.31657 555.31988
61 154.78646 -621.31657
62 2164.10178 154.78646
63 736.89615 2164.10178
64 484.28046 736.89615
65 2355.06993 484.28046
66 -1793.04204 2355.06993
67 -1240.46086 -1793.04204
68 1472.59406 -1240.46086
69 266.20796 1472.59406
70 2293.73589 266.20796
71 2339.85884 2293.73589
72 1975.69649 2339.85884
73 1797.90931 1975.69649
74 5289.67732 1797.90931
75 1200.62482 5289.67732
76 2937.87665 1200.62482
77 3124.18885 2937.87665
78 69.71715 3124.18885
79 559.25992 69.71715
80 2466.84473 559.25992
81 3613.49623 2466.84473
82 3855.54559 3613.49623
83 2247.27201 3855.54559
> plot(z,main=paste('Residual Lag plot, lowess, and regression line'), ylab='values of Residuals', xlab='lagged values of Residuals')
> lines(lowess(z))
> abline(lm(z))
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/7ay3n1229261596.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> acf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Autocorrelation Function')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/8f0w01229261596.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> pacf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Partial Autocorrelation Function')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/9s6861229261596.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> opar <- par(mfrow = c(2,2), oma = c(0, 0, 1.1, 0))
> plot(mylm, las = 1, sub='Residual Diagnostics')
> par(opar)
> dev.off()
null device
1
> if (n > n25) {
+ postscript(file="/var/www/html/freestat/rcomp/tmp/100jnf1229261596.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
+ plot(kp3:nmkm3,gqarr[,2], main='Goldfeld-Quandt test',ylab='2-sided p-value',xlab='breakpoint')
+ grid()
+ dev.off()
+ }
null device
1
>
> #Note: the /var/www/html/freestat/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/html/freestat/rcomp/createtable")
>
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Estimated Regression Equation', 1, TRUE)
> a<-table.row.end(a)
> myeq <- colnames(x)[1]
> myeq <- paste(myeq, '[t] = ', sep='')
> for (i in 1:k){
+ if (mysum$coefficients[i,1] > 0) myeq <- paste(myeq, '+', '')
+ myeq <- paste(myeq, mysum$coefficients[i,1], sep=' ')
+ if (rownames(mysum$coefficients)[i] != '(Intercept)') {
+ myeq <- paste(myeq, rownames(mysum$coefficients)[i], sep='')
+ if (rownames(mysum$coefficients)[i] != 't') myeq <- paste(myeq, '[t]', sep='')
+ }
+ }
> myeq <- paste(myeq, ' + e[t]')
> a<-table.row.start(a)
> a<-table.element(a, myeq)
> a<-table.row.end(a)
> a<-table.end(a)
> table.save(a,file="/var/www/html/freestat/rcomp/tmp/1125xb1229261596.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a,hyperlink('http://www.xycoon.com/ols1.htm','Multiple Linear Regression - Ordinary Least Squares',''), 6, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a,'Variable',header=TRUE)
> a<-table.element(a,'Parameter',header=TRUE)
> a<-table.element(a,'S.D.',header=TRUE)
> a<-table.element(a,'T-STAT
H0: parameter = 0',header=TRUE)
> a<-table.element(a,'2-tail p-value',header=TRUE)
> a<-table.element(a,'1-tail p-value',header=TRUE)
> a<-table.row.end(a)
> for (i in 1:k){
+ a<-table.row.start(a)
+ a<-table.element(a,rownames(mysum$coefficients)[i],header=TRUE)
+ a<-table.element(a,mysum$coefficients[i,1])
+ a<-table.element(a, round(mysum$coefficients[i,2],6))
+ a<-table.element(a, round(mysum$coefficients[i,3],4))
+ a<-table.element(a, round(mysum$coefficients[i,4],6))
+ a<-table.element(a, round(mysum$coefficients[i,4]/2,6))
+ a<-table.row.end(a)
+ }
> a<-table.end(a)
> table.save(a,file="/var/www/html/freestat/rcomp/tmp/12nb401229261596.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Regression Statistics', 2, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple R',1,TRUE)
> a<-table.element(a, sqrt(mysum$r.squared))
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'R-squared',1,TRUE)
> a<-table.element(a, mysum$r.squared)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Adjusted R-squared',1,TRUE)
> a<-table.element(a, mysum$adj.r.squared)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (value)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[1])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (DF numerator)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[2])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (DF denominator)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[3])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'p-value',1,TRUE)
> a<-table.element(a, 1-pf(mysum$fstatistic[1],mysum$fstatistic[2],mysum$fstatistic[3]))
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Residual Statistics', 2, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Residual Standard Deviation',1,TRUE)
> a<-table.element(a, mysum$sigma)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Sum Squared Residuals',1,TRUE)
> a<-table.element(a, sum(myerror*myerror))
> a<-table.row.end(a)
> a<-table.end(a)
> table.save(a,file="/var/www/html/freestat/rcomp/tmp/133nos1229261597.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Actuals, Interpolation, and Residuals', 4, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Time or Index', 1, TRUE)
> a<-table.element(a, 'Actuals', 1, TRUE)
> a<-table.element(a, 'Interpolation
Forecast', 1, TRUE)
> a<-table.element(a, 'Residuals
Prediction Error', 1, TRUE)
> a<-table.row.end(a)
> for (i in 1:n) {
+ a<-table.row.start(a)
+ a<-table.element(a,i, 1, TRUE)
+ a<-table.element(a,x[i])
+ a<-table.element(a,x[i]-mysum$resid[i])
+ a<-table.element(a,mysum$resid[i])
+ a<-table.row.end(a)
+ }
> a<-table.end(a)
> table.save(a,file="/var/www/html/freestat/rcomp/tmp/14lt0e1229261597.tab")
> if (n > n25) {
+ a<-table.start()
+ a<-table.row.start(a)
+ a<-table.element(a,'Goldfeld-Quandt test for Heteroskedasticity',4,TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'p-values',header=TRUE)
+ a<-table.element(a,'Alternative Hypothesis',3,header=TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'breakpoint index',header=TRUE)
+ a<-table.element(a,'greater',header=TRUE)
+ a<-table.element(a,'2-sided',header=TRUE)
+ a<-table.element(a,'less',header=TRUE)
+ a<-table.row.end(a)
+ for (mypoint in kp3:nmkm3) {
+ a<-table.row.start(a)
+ a<-table.element(a,mypoint,header=TRUE)
+ a<-table.element(a,gqarr[mypoint-kp3+1,1])
+ a<-table.element(a,gqarr[mypoint-kp3+1,2])
+ a<-table.element(a,gqarr[mypoint-kp3+1,3])
+ a<-table.row.end(a)
+ }
+ a<-table.end(a)
+ table.save(a,file="/var/www/html/freestat/rcomp/tmp/15hr091229261597.tab")
+ a<-table.start()
+ a<-table.row.start(a)
+ a<-table.element(a,'Meta Analysis of Goldfeld-Quandt test for Heteroskedasticity',4,TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'Description',header=TRUE)
+ a<-table.element(a,'# significant tests',header=TRUE)
+ a<-table.element(a,'% significant tests',header=TRUE)
+ a<-table.element(a,'OK/NOK',header=TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'1% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant1)
+ a<-table.element(a,numsignificant1/numgqtests)
+ if (numsignificant1/numgqtests < 0.01) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'5% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant5)
+ a<-table.element(a,numsignificant5/numgqtests)
+ if (numsignificant5/numgqtests < 0.05) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'10% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant10)
+ a<-table.element(a,numsignificant10/numgqtests)
+ if (numsignificant10/numgqtests < 0.1) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.end(a)
+ table.save(a,file="/var/www/html/freestat/rcomp/tmp/16uez81229261597.tab")
+ }
>
> system("convert tmp/1sauh1229261596.ps tmp/1sauh1229261596.png")
> system("convert tmp/2o1hb1229261596.ps tmp/2o1hb1229261596.png")
> system("convert tmp/3edx11229261596.ps tmp/3edx11229261596.png")
> system("convert tmp/4h1d51229261596.ps tmp/4h1d51229261596.png")
> system("convert tmp/5nno81229261596.ps tmp/5nno81229261596.png")
> system("convert tmp/65cgk1229261596.ps tmp/65cgk1229261596.png")
> system("convert tmp/7ay3n1229261596.ps tmp/7ay3n1229261596.png")
> system("convert tmp/8f0w01229261596.ps tmp/8f0w01229261596.png")
> system("convert tmp/9s6861229261596.ps tmp/9s6861229261596.png")
> system("convert tmp/100jnf1229261596.ps tmp/100jnf1229261596.png")
>
>
> proc.time()
user system elapsed
3.964 2.490 4.336