R version 2.9.0 (2009-04-17)
Copyright (C) 2009 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(100
+ ,0
+ ,97.56592292
+ ,0
+ ,93.71196755
+ ,0
+ ,92.69776876
+ ,0
+ ,89.65517241
+ ,0
+ ,89.04665314
+ ,0
+ ,98.98580122
+ ,0
+ ,105.6795132
+ ,0
+ ,101.6227181
+ ,0
+ ,98.37728195
+ ,0
+ ,94.11764706
+ ,0
+ ,93.30628803
+ ,0
+ ,94.72616633
+ ,0
+ ,93.30628803
+ ,0
+ ,90.87221095
+ ,0
+ ,89.85801217
+ ,0
+ ,88.43813387
+ ,0
+ ,87.42393509
+ ,0
+ ,98.17444219
+ ,0
+ ,103.4482759
+ ,0
+ ,104.0567951
+ ,0
+ ,102.0283976
+ ,0
+ ,95.53752535
+ ,0
+ ,95.53752535
+ ,0
+ ,96.55172414
+ ,0
+ ,96.34888438
+ ,0
+ ,95.3346856
+ ,0
+ ,93.50912779
+ ,0
+ ,92.29208925
+ ,0
+ ,92.49492901
+ ,0
+ ,104.8681542
+ ,0
+ ,106.4908722
+ ,0
+ ,106.0851927
+ ,0
+ ,105.2738337
+ ,0
+ ,103.2454361
+ ,0
+ ,103.8539554
+ ,0
+ ,105.2738337
+ ,0
+ ,104.8681542
+ ,0
+ ,103.4482759
+ ,0
+ ,103.2454361
+ ,0
+ ,101.6227181
+ ,0
+ ,102.8397566
+ ,0
+ ,115.4158215
+ ,0
+ ,117.6470588
+ ,0
+ ,117.2413793
+ ,0
+ ,114.6044625
+ ,0
+ ,110.9533469
+ ,0
+ ,112.5760649
+ ,0
+ ,113.9959432
+ ,0
+ ,113.7931034
+ ,0
+ ,112.5760649
+ ,0
+ ,110.3448276
+ ,0
+ ,108.9249493
+ ,0
+ ,110.1419878
+ ,0
+ ,120.4868154
+ ,0
+ ,123.9350913
+ ,0
+ ,124.3407708
+ ,0
+ ,123.9350913
+ ,0
+ ,120.4868154
+ ,0
+ ,120.6896552
+ ,0
+ ,119.8782961
+ ,0
+ ,119.4726166
+ ,0
+ ,118.4584178
+ ,0
+ ,116.2271805
+ ,0
+ ,115.010142
+ ,0
+ ,115.4158215
+ ,0
+ ,125.9634888
+ ,0
+ ,127.5862069
+ ,0
+ ,127.3833671
+ ,0
+ ,124.137931
+ ,0
+ ,120.6896552
+ ,0
+ ,121.0953347
+ ,0
+ ,120.2839757
+ ,0
+ ,119.6754564
+ ,0
+ ,117.6470588
+ ,0
+ ,116.4300203
+ ,0
+ ,116.2271805
+ ,0
+ ,116.2271805
+ ,0
+ ,125.7606491
+ ,0
+ ,126.9776876
+ ,0
+ ,125.7606491
+ ,0
+ ,119.2697769
+ ,0
+ ,114.8073022
+ ,0
+ ,112.9817444
+ ,0
+ ,113.7931034
+ ,0
+ ,111.3590264
+ ,0
+ ,107.9107505
+ ,0
+ ,106.693712
+ ,0
+ ,103.6511156
+ ,0
+ ,101.2170385
+ ,0
+ ,112.5760649
+ ,0
+ ,114.6044625
+ ,1
+ ,109.9391481
+ ,1
+ ,106.8965517
+ ,1
+ ,103.4482759
+ ,1
+ ,104.2596349
+ ,1
+ ,104.8681542
+ ,1
+ ,103.0425963
+ ,1
+ ,100
+ ,1
+ ,99.39148073
+ ,1
+ ,95.13184584
+ ,1
+ ,96.95740365
+ ,1
+ ,107.0993915
+ ,1
+ ,108.31643
+ ,1
+ ,105.0709939
+ ,1
+ ,102.6369168
+ ,1
+ ,101.8255578
+ ,1
+ ,104.6653144
+ ,1)
+ ,dim=c(2
+ ,108)
+ ,dimnames=list(c('TW'
+ ,'DUM')
+ ,1:108))
> y <- array(NA,dim=c(2,108),dimnames=list(c('TW','DUM'),1:108))
> 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
TW DUM
1 100.00000 0
2 97.56592 0
3 93.71197 0
4 92.69777 0
5 89.65517 0
6 89.04665 0
7 98.98580 0
8 105.67951 0
9 101.62272 0
10 98.37728 0
11 94.11765 0
12 93.30629 0
13 94.72617 0
14 93.30629 0
15 90.87221 0
16 89.85801 0
17 88.43813 0
18 87.42394 0
19 98.17444 0
20 103.44828 0
21 104.05680 0
22 102.02840 0
23 95.53753 0
24 95.53753 0
25 96.55172 0
26 96.34888 0
27 95.33469 0
28 93.50913 0
29 92.29209 0
30 92.49493 0
31 104.86815 0
32 106.49087 0
33 106.08519 0
34 105.27383 0
35 103.24544 0
36 103.85396 0
37 105.27383 0
38 104.86815 0
39 103.44828 0
40 103.24544 0
41 101.62272 0
42 102.83976 0
43 115.41582 0
44 117.64706 0
45 117.24138 0
46 114.60446 0
47 110.95335 0
48 112.57606 0
49 113.99594 0
50 113.79310 0
51 112.57606 0
52 110.34483 0
53 108.92495 0
54 110.14199 0
55 120.48682 0
56 123.93509 0
57 124.34077 0
58 123.93509 0
59 120.48682 0
60 120.68966 0
61 119.87830 0
62 119.47262 0
63 118.45842 0
64 116.22718 0
65 115.01014 0
66 115.41582 0
67 125.96349 0
68 127.58621 0
69 127.38337 0
70 124.13793 0
71 120.68966 0
72 121.09533 0
73 120.28398 0
74 119.67546 0
75 117.64706 0
76 116.43002 0
77 116.22718 0
78 116.22718 0
79 125.76065 0
80 126.97769 0
81 125.76065 0
82 119.26978 0
83 114.80730 0
84 112.98174 0
85 113.79310 0
86 111.35903 0
87 107.91075 0
88 106.69371 0
89 103.65112 0
90 101.21704 0
91 112.57606 0
92 114.60446 1
93 109.93915 1
94 106.89655 1
95 103.44828 1
96 104.25963 1
97 104.86815 1
98 103.04260 1
99 100.00000 1
100 99.39148 1
101 95.13185 1
102 96.95740 1
103 107.09939 1
104 108.31643 1
105 105.07099 1
106 102.63692 1
107 101.82556 1
108 104.66531 1
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) DUM
108.225 -4.216
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-20.80111 -7.01892 -0.03186 8.00214 19.36117
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 108.225 1.101 98.272 <2e-16 ***
DUM -4.216 2.776 -1.519 0.132
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 10.51 on 106 degrees of freedom
Multiple R-squared: 0.0213, Adjusted R-squared: 0.01207
F-statistic: 2.307 on 1 and 106 DF, p-value: 0.1318
> 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.111150607 2.223012e-01 8.888494e-01
[2,] 0.077875465 1.557509e-01 9.221245e-01
[3,] 0.050352442 1.007049e-01 9.496476e-01
[4,] 0.099639432 1.992789e-01 9.003606e-01
[5,] 0.071504368 1.430087e-01 9.284956e-01
[6,] 0.039608236 7.921647e-02 9.603918e-01
[7,] 0.023099272 4.619854e-02 9.769007e-01
[8,] 0.014168267 2.833653e-02 9.858317e-01
[9,] 0.007704948 1.540990e-02 9.922951e-01
[10,] 0.004652694 9.305389e-03 9.953473e-01
[11,] 0.003973833 7.947666e-03 9.960262e-01
[12,] 0.004069324 8.138648e-03 9.959307e-01
[13,] 0.005492606 1.098521e-02 9.945074e-01
[14,] 0.008911873 1.782375e-02 9.910881e-01
[15,] 0.006711382 1.342276e-02 9.932886e-01
[16,] 0.009685302 1.937060e-02 9.903147e-01
[17,] 0.013314816 2.662963e-02 9.866852e-01
[18,] 0.012907538 2.581508e-02 9.870925e-01
[19,] 0.010149762 2.029952e-02 9.898502e-01
[20,] 0.008225204 1.645041e-02 9.917748e-01
[21,] 0.006676639 1.335328e-02 9.933234e-01
[22,] 0.005647649 1.129530e-02 9.943524e-01
[23,] 0.005250568 1.050114e-02 9.947494e-01
[24,] 0.006144084 1.228817e-02 9.938559e-01
[25,] 0.009307353 1.861471e-02 9.906926e-01
[26,] 0.015423304 3.084661e-02 9.845767e-01
[27,] 0.027917980 5.583596e-02 9.720820e-01
[28,] 0.053004773 1.060095e-01 9.469952e-01
[29,] 0.082121571 1.642431e-01 9.178784e-01
[30,] 0.109824621 2.196492e-01 8.901754e-01
[31,] 0.130618188 2.612364e-01 8.693818e-01
[32,] 0.157274489 3.145490e-01 8.427255e-01
[33,] 0.194809455 3.896189e-01 8.051905e-01
[34,] 0.233974355 4.679487e-01 7.660256e-01
[35,] 0.275443846 5.508877e-01 7.245562e-01
[36,] 0.326670748 6.533415e-01 6.733293e-01
[37,] 0.401041036 8.020821e-01 5.989590e-01
[38,] 0.486093279 9.721866e-01 5.139067e-01
[39,] 0.700221851 5.995563e-01 2.997781e-01
[40,] 0.862519112 2.749618e-01 1.374809e-01
[41,] 0.932858469 1.342831e-01 6.714153e-02
[42,] 0.955798771 8.840246e-02 4.420123e-02
[43,] 0.964001478 7.199704e-02 3.599852e-02
[44,] 0.971360671 5.727866e-02 2.863933e-02
[45,] 0.977764647 4.447071e-02 2.223535e-02
[46,] 0.981814285 3.637143e-02 1.818572e-02
[47,] 0.984032485 3.193503e-02 1.596752e-02
[48,] 0.985725501 2.854900e-02 1.427450e-02
[49,] 0.987905516 2.418897e-02 1.209448e-02
[50,] 0.989498711 2.100258e-02 1.050129e-02
[51,] 0.993937675 1.212465e-02 6.062325e-03
[52,] 0.997665948 4.668103e-03 2.334052e-03
[53,] 0.999083517 1.832967e-03 9.164833e-04
[54,] 0.999578189 8.436223e-04 4.218111e-04
[55,] 0.999652137 6.957266e-04 3.478633e-04
[56,] 0.999702529 5.949427e-04 2.974714e-04
[57,] 0.999708184 5.836319e-04 2.918159e-04
[58,] 0.999690366 6.192671e-04 3.096336e-04
[59,] 0.999633345 7.333098e-04 3.666549e-04
[60,] 0.999509484 9.810325e-04 4.905162e-04
[61,] 0.999327331 1.345338e-03 6.726690e-04
[62,] 0.999077483 1.845034e-03 9.225168e-04
[63,] 0.999496367 1.007265e-03 5.036326e-04
[64,] 0.999805068 3.898649e-04 1.949324e-04
[65,] 0.999927844 1.443119e-04 7.215597e-05
[66,] 0.999950151 9.969818e-05 4.984909e-05
[67,] 0.999941288 1.174235e-04 5.871174e-05
[68,] 0.999934653 1.306935e-04 6.534673e-05
[69,] 0.999919859 1.602828e-04 8.014141e-05
[70,] 0.999895400 2.092008e-04 1.046004e-04
[71,] 0.999836154 3.276924e-04 1.638462e-04
[72,] 0.999725461 5.490789e-04 2.745395e-04
[73,] 0.999540882 9.182354e-04 4.591177e-04
[74,] 0.999242483 1.515035e-03 7.575175e-04
[75,] 0.999677946 6.441087e-04 3.220543e-04
[76,] 0.999941183 1.176341e-04 5.881703e-05
[77,] 0.999994127 1.174608e-05 5.873041e-06
[78,] 0.999996952 6.096907e-06 3.048454e-06
[79,] 0.999995809 8.382403e-06 4.191202e-06
[80,] 0.999992702 1.459596e-05 7.297980e-06
[81,] 0.999990766 1.846840e-05 9.234200e-06
[82,] 0.999983898 3.220413e-05 1.610206e-05
[83,] 0.999961053 7.789450e-05 3.894725e-05
[84,] 0.999904619 1.907628e-04 9.538139e-05
[85,] 0.999797680 4.046393e-04 2.023196e-04
[86,] 0.999841252 3.174959e-04 1.587479e-04
[87,] 0.999612716 7.745672e-04 3.872836e-04
[88,] 0.999912805 1.743909e-04 8.719544e-05
[89,] 0.999927968 1.440643e-04 7.203213e-05
[90,] 0.999869172 2.616567e-04 1.308284e-04
[91,] 0.999622766 7.544677e-04 3.772338e-04
[92,] 0.999005960 1.988080e-03 9.940398e-04
[93,] 0.997647542 4.704916e-03 2.352458e-03
[94,] 0.993862307 1.227539e-02 6.137693e-03
[95,] 0.985976047 2.804791e-02 1.402395e-02
[96,] 0.971366309 5.726738e-02 2.863369e-02
[97,] 0.980756487 3.848703e-02 1.924351e-02
[98,] 0.991694752 1.661050e-02 8.305248e-03
[99,] 0.973983105 5.203379e-02 2.601689e-02
> postscript(file="/var/www/html/rcomp/tmp/1nc401260803020.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/rcomp/tmp/267a71260803020.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/rcomp/tmp/3p76w1260803020.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/rcomp/tmp/4bmzu1260803020.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/rcomp/tmp/5dkcm1260803020.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 = 108
Frequency = 1
1 2 3 4 5 6
-8.2250407 -10.6591178 -14.5130731 -15.5272719 -18.5698683 -19.1783875
7 8 9 10 11 12
-9.2392395 -2.5455275 -6.6023226 -9.8477587 -14.1073936 -14.9187526
13 14 15 16 17 18
-13.4988743 -14.9187526 -17.3528297 -18.3670285 -19.7869068 -20.8011056
19 20 21 22 23 24
-10.0505985 -4.7767648 -4.1682456 -6.1966431 -12.6875153 -12.6875153
25 26 27 28 29 30
-11.6733165 -11.8761563 -12.8903551 -14.7159129 -15.9329514 -15.7301117
31 32 33 34 35 36
-3.3568865 -1.7341685 -2.1398480 -2.9512070 -4.9796046 -4.3710853
37 38 39 40 41 42
-2.9512070 -3.3568865 -4.7767648 -4.9796046 -6.6023226 -5.3852841
43 44 45 46 47 48
7.1907808 9.4220181 9.0163386 6.3794218 2.7283062 4.3510242
49 50 51 52 53 54
5.7709025 5.5680627 4.3510242 2.1197869 0.6999086 1.9169471
55 56 57 58 59 60
12.2617747 15.7100506 16.1157301 15.7100506 12.2617747 12.4646145
61 62 63 64 65 66
11.6532554 11.2475759 10.2333771 8.0021398 6.7851013 7.1907808
67 68 69 70 71 72
17.7384481 19.3611662 19.1583264 15.9128903 12.4646145 12.8702940
73 74 75 76 77 78
12.0589350 11.4504157 9.4220181 8.2049796 8.0021398 8.0021398
79 80 81 82 83 84
17.5356084 18.7526469 17.5356084 11.0447362 6.5822615 4.7567037
85 86 87 88 89 90
5.5680627 3.1339857 -0.3142902 -1.5313287 -4.5739251 -7.0080022
91 92 93 94 95 96
4.3510242 10.5953944 5.9300800 2.8874836 -0.5607922 0.2505668
97 98 99 100 101 102
0.8590861 -0.9664718 -4.0090681 -4.6175874 -8.8772223 -7.0516645
103 104 105 106 107 108
3.0903234 4.3073619 1.0619258 -1.3721513 -2.1835103 0.6562463
> postscript(file="/var/www/html/rcomp/tmp/6luah1260803020.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 = 108
Frequency = 1
lag(myerror, k = 1) myerror
0 -8.2250407 NA
1 -10.6591178 -8.2250407
2 -14.5130731 -10.6591178
3 -15.5272719 -14.5130731
4 -18.5698683 -15.5272719
5 -19.1783875 -18.5698683
6 -9.2392395 -19.1783875
7 -2.5455275 -9.2392395
8 -6.6023226 -2.5455275
9 -9.8477587 -6.6023226
10 -14.1073936 -9.8477587
11 -14.9187526 -14.1073936
12 -13.4988743 -14.9187526
13 -14.9187526 -13.4988743
14 -17.3528297 -14.9187526
15 -18.3670285 -17.3528297
16 -19.7869068 -18.3670285
17 -20.8011056 -19.7869068
18 -10.0505985 -20.8011056
19 -4.7767648 -10.0505985
20 -4.1682456 -4.7767648
21 -6.1966431 -4.1682456
22 -12.6875153 -6.1966431
23 -12.6875153 -12.6875153
24 -11.6733165 -12.6875153
25 -11.8761563 -11.6733165
26 -12.8903551 -11.8761563
27 -14.7159129 -12.8903551
28 -15.9329514 -14.7159129
29 -15.7301117 -15.9329514
30 -3.3568865 -15.7301117
31 -1.7341685 -3.3568865
32 -2.1398480 -1.7341685
33 -2.9512070 -2.1398480
34 -4.9796046 -2.9512070
35 -4.3710853 -4.9796046
36 -2.9512070 -4.3710853
37 -3.3568865 -2.9512070
38 -4.7767648 -3.3568865
39 -4.9796046 -4.7767648
40 -6.6023226 -4.9796046
41 -5.3852841 -6.6023226
42 7.1907808 -5.3852841
43 9.4220181 7.1907808
44 9.0163386 9.4220181
45 6.3794218 9.0163386
46 2.7283062 6.3794218
47 4.3510242 2.7283062
48 5.7709025 4.3510242
49 5.5680627 5.7709025
50 4.3510242 5.5680627
51 2.1197869 4.3510242
52 0.6999086 2.1197869
53 1.9169471 0.6999086
54 12.2617747 1.9169471
55 15.7100506 12.2617747
56 16.1157301 15.7100506
57 15.7100506 16.1157301
58 12.2617747 15.7100506
59 12.4646145 12.2617747
60 11.6532554 12.4646145
61 11.2475759 11.6532554
62 10.2333771 11.2475759
63 8.0021398 10.2333771
64 6.7851013 8.0021398
65 7.1907808 6.7851013
66 17.7384481 7.1907808
67 19.3611662 17.7384481
68 19.1583264 19.3611662
69 15.9128903 19.1583264
70 12.4646145 15.9128903
71 12.8702940 12.4646145
72 12.0589350 12.8702940
73 11.4504157 12.0589350
74 9.4220181 11.4504157
75 8.2049796 9.4220181
76 8.0021398 8.2049796
77 8.0021398 8.0021398
78 17.5356084 8.0021398
79 18.7526469 17.5356084
80 17.5356084 18.7526469
81 11.0447362 17.5356084
82 6.5822615 11.0447362
83 4.7567037 6.5822615
84 5.5680627 4.7567037
85 3.1339857 5.5680627
86 -0.3142902 3.1339857
87 -1.5313287 -0.3142902
88 -4.5739251 -1.5313287
89 -7.0080022 -4.5739251
90 4.3510242 -7.0080022
91 10.5953944 4.3510242
92 5.9300800 10.5953944
93 2.8874836 5.9300800
94 -0.5607922 2.8874836
95 0.2505668 -0.5607922
96 0.8590861 0.2505668
97 -0.9664718 0.8590861
98 -4.0090681 -0.9664718
99 -4.6175874 -4.0090681
100 -8.8772223 -4.6175874
101 -7.0516645 -8.8772223
102 3.0903234 -7.0516645
103 4.3073619 3.0903234
104 1.0619258 4.3073619
105 -1.3721513 1.0619258
106 -2.1835103 -1.3721513
107 0.6562463 -2.1835103
108 NA 0.6562463
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -10.6591178 -8.2250407
[2,] -14.5130731 -10.6591178
[3,] -15.5272719 -14.5130731
[4,] -18.5698683 -15.5272719
[5,] -19.1783875 -18.5698683
[6,] -9.2392395 -19.1783875
[7,] -2.5455275 -9.2392395
[8,] -6.6023226 -2.5455275
[9,] -9.8477587 -6.6023226
[10,] -14.1073936 -9.8477587
[11,] -14.9187526 -14.1073936
[12,] -13.4988743 -14.9187526
[13,] -14.9187526 -13.4988743
[14,] -17.3528297 -14.9187526
[15,] -18.3670285 -17.3528297
[16,] -19.7869068 -18.3670285
[17,] -20.8011056 -19.7869068
[18,] -10.0505985 -20.8011056
[19,] -4.7767648 -10.0505985
[20,] -4.1682456 -4.7767648
[21,] -6.1966431 -4.1682456
[22,] -12.6875153 -6.1966431
[23,] -12.6875153 -12.6875153
[24,] -11.6733165 -12.6875153
[25,] -11.8761563 -11.6733165
[26,] -12.8903551 -11.8761563
[27,] -14.7159129 -12.8903551
[28,] -15.9329514 -14.7159129
[29,] -15.7301117 -15.9329514
[30,] -3.3568865 -15.7301117
[31,] -1.7341685 -3.3568865
[32,] -2.1398480 -1.7341685
[33,] -2.9512070 -2.1398480
[34,] -4.9796046 -2.9512070
[35,] -4.3710853 -4.9796046
[36,] -2.9512070 -4.3710853
[37,] -3.3568865 -2.9512070
[38,] -4.7767648 -3.3568865
[39,] -4.9796046 -4.7767648
[40,] -6.6023226 -4.9796046
[41,] -5.3852841 -6.6023226
[42,] 7.1907808 -5.3852841
[43,] 9.4220181 7.1907808
[44,] 9.0163386 9.4220181
[45,] 6.3794218 9.0163386
[46,] 2.7283062 6.3794218
[47,] 4.3510242 2.7283062
[48,] 5.7709025 4.3510242
[49,] 5.5680627 5.7709025
[50,] 4.3510242 5.5680627
[51,] 2.1197869 4.3510242
[52,] 0.6999086 2.1197869
[53,] 1.9169471 0.6999086
[54,] 12.2617747 1.9169471
[55,] 15.7100506 12.2617747
[56,] 16.1157301 15.7100506
[57,] 15.7100506 16.1157301
[58,] 12.2617747 15.7100506
[59,] 12.4646145 12.2617747
[60,] 11.6532554 12.4646145
[61,] 11.2475759 11.6532554
[62,] 10.2333771 11.2475759
[63,] 8.0021398 10.2333771
[64,] 6.7851013 8.0021398
[65,] 7.1907808 6.7851013
[66,] 17.7384481 7.1907808
[67,] 19.3611662 17.7384481
[68,] 19.1583264 19.3611662
[69,] 15.9128903 19.1583264
[70,] 12.4646145 15.9128903
[71,] 12.8702940 12.4646145
[72,] 12.0589350 12.8702940
[73,] 11.4504157 12.0589350
[74,] 9.4220181 11.4504157
[75,] 8.2049796 9.4220181
[76,] 8.0021398 8.2049796
[77,] 8.0021398 8.0021398
[78,] 17.5356084 8.0021398
[79,] 18.7526469 17.5356084
[80,] 17.5356084 18.7526469
[81,] 11.0447362 17.5356084
[82,] 6.5822615 11.0447362
[83,] 4.7567037 6.5822615
[84,] 5.5680627 4.7567037
[85,] 3.1339857 5.5680627
[86,] -0.3142902 3.1339857
[87,] -1.5313287 -0.3142902
[88,] -4.5739251 -1.5313287
[89,] -7.0080022 -4.5739251
[90,] 4.3510242 -7.0080022
[91,] 10.5953944 4.3510242
[92,] 5.9300800 10.5953944
[93,] 2.8874836 5.9300800
[94,] -0.5607922 2.8874836
[95,] 0.2505668 -0.5607922
[96,] 0.8590861 0.2505668
[97,] -0.9664718 0.8590861
[98,] -4.0090681 -0.9664718
[99,] -4.6175874 -4.0090681
[100,] -8.8772223 -4.6175874
[101,] -7.0516645 -8.8772223
[102,] 3.0903234 -7.0516645
[103,] 4.3073619 3.0903234
[104,] 1.0619258 4.3073619
[105,] -1.3721513 1.0619258
[106,] -2.1835103 -1.3721513
[107,] 0.6562463 -2.1835103
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -10.6591178 -8.2250407
2 -14.5130731 -10.6591178
3 -15.5272719 -14.5130731
4 -18.5698683 -15.5272719
5 -19.1783875 -18.5698683
6 -9.2392395 -19.1783875
7 -2.5455275 -9.2392395
8 -6.6023226 -2.5455275
9 -9.8477587 -6.6023226
10 -14.1073936 -9.8477587
11 -14.9187526 -14.1073936
12 -13.4988743 -14.9187526
13 -14.9187526 -13.4988743
14 -17.3528297 -14.9187526
15 -18.3670285 -17.3528297
16 -19.7869068 -18.3670285
17 -20.8011056 -19.7869068
18 -10.0505985 -20.8011056
19 -4.7767648 -10.0505985
20 -4.1682456 -4.7767648
21 -6.1966431 -4.1682456
22 -12.6875153 -6.1966431
23 -12.6875153 -12.6875153
24 -11.6733165 -12.6875153
25 -11.8761563 -11.6733165
26 -12.8903551 -11.8761563
27 -14.7159129 -12.8903551
28 -15.9329514 -14.7159129
29 -15.7301117 -15.9329514
30 -3.3568865 -15.7301117
31 -1.7341685 -3.3568865
32 -2.1398480 -1.7341685
33 -2.9512070 -2.1398480
34 -4.9796046 -2.9512070
35 -4.3710853 -4.9796046
36 -2.9512070 -4.3710853
37 -3.3568865 -2.9512070
38 -4.7767648 -3.3568865
39 -4.9796046 -4.7767648
40 -6.6023226 -4.9796046
41 -5.3852841 -6.6023226
42 7.1907808 -5.3852841
43 9.4220181 7.1907808
44 9.0163386 9.4220181
45 6.3794218 9.0163386
46 2.7283062 6.3794218
47 4.3510242 2.7283062
48 5.7709025 4.3510242
49 5.5680627 5.7709025
50 4.3510242 5.5680627
51 2.1197869 4.3510242
52 0.6999086 2.1197869
53 1.9169471 0.6999086
54 12.2617747 1.9169471
55 15.7100506 12.2617747
56 16.1157301 15.7100506
57 15.7100506 16.1157301
58 12.2617747 15.7100506
59 12.4646145 12.2617747
60 11.6532554 12.4646145
61 11.2475759 11.6532554
62 10.2333771 11.2475759
63 8.0021398 10.2333771
64 6.7851013 8.0021398
65 7.1907808 6.7851013
66 17.7384481 7.1907808
67 19.3611662 17.7384481
68 19.1583264 19.3611662
69 15.9128903 19.1583264
70 12.4646145 15.9128903
71 12.8702940 12.4646145
72 12.0589350 12.8702940
73 11.4504157 12.0589350
74 9.4220181 11.4504157
75 8.2049796 9.4220181
76 8.0021398 8.2049796
77 8.0021398 8.0021398
78 17.5356084 8.0021398
79 18.7526469 17.5356084
80 17.5356084 18.7526469
81 11.0447362 17.5356084
82 6.5822615 11.0447362
83 4.7567037 6.5822615
84 5.5680627 4.7567037
85 3.1339857 5.5680627
86 -0.3142902 3.1339857
87 -1.5313287 -0.3142902
88 -4.5739251 -1.5313287
89 -7.0080022 -4.5739251
90 4.3510242 -7.0080022
91 10.5953944 4.3510242
92 5.9300800 10.5953944
93 2.8874836 5.9300800
94 -0.5607922 2.8874836
95 0.2505668 -0.5607922
96 0.8590861 0.2505668
97 -0.9664718 0.8590861
98 -4.0090681 -0.9664718
99 -4.6175874 -4.0090681
100 -8.8772223 -4.6175874
101 -7.0516645 -8.8772223
102 3.0903234 -7.0516645
103 4.3073619 3.0903234
104 1.0619258 4.3073619
105 -1.3721513 1.0619258
106 -2.1835103 -1.3721513
107 0.6562463 -2.1835103
> 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/rcomp/tmp/755cp1260803020.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/rcomp/tmp/8inh81260803020.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/rcomp/tmp/9wtr31260803020.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/rcomp/tmp/10jieb1260803020.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/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/html/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/rcomp/tmp/11s7cj1260803020.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/rcomp/tmp/125atz1260803020.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/rcomp/tmp/133vve1260803020.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/rcomp/tmp/14fwer1260803020.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/rcomp/tmp/15yj9l1260803020.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/rcomp/tmp/164yvg1260803020.tab")
+ }
>
> try(system("convert tmp/1nc401260803020.ps tmp/1nc401260803020.png",intern=TRUE))
character(0)
> try(system("convert tmp/267a71260803020.ps tmp/267a71260803020.png",intern=TRUE))
character(0)
> try(system("convert tmp/3p76w1260803020.ps tmp/3p76w1260803020.png",intern=TRUE))
character(0)
> try(system("convert tmp/4bmzu1260803020.ps tmp/4bmzu1260803020.png",intern=TRUE))
character(0)
> try(system("convert tmp/5dkcm1260803020.ps tmp/5dkcm1260803020.png",intern=TRUE))
character(0)
> try(system("convert tmp/6luah1260803020.ps tmp/6luah1260803020.png",intern=TRUE))
character(0)
> try(system("convert tmp/755cp1260803020.ps tmp/755cp1260803020.png",intern=TRUE))
character(0)
> try(system("convert tmp/8inh81260803020.ps tmp/8inh81260803020.png",intern=TRUE))
character(0)
> try(system("convert tmp/9wtr31260803020.ps tmp/9wtr31260803020.png",intern=TRUE))
character(0)
> try(system("convert tmp/10jieb1260803020.ps tmp/10jieb1260803020.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
2.980 1.587 4.502