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(4716.99
+ ,0
+ ,4926.65
+ ,0
+ ,4920.10
+ ,0
+ ,5170.09
+ ,0
+ ,5246.24
+ ,0
+ ,5283.61
+ ,0
+ ,4979.05
+ ,0
+ ,4825.20
+ ,0
+ ,4695.12
+ ,0
+ ,4711.54
+ ,0
+ ,4727.22
+ ,0
+ ,4384.96
+ ,0
+ ,4378.75
+ ,0
+ ,4472.93
+ ,0
+ ,4564.07
+ ,0
+ ,4310.54
+ ,0
+ ,4171.38
+ ,0
+ ,4049.38
+ ,0
+ ,3591.37
+ ,0
+ ,3720.46
+ ,0
+ ,4107.23
+ ,0
+ ,4101.71
+ ,0
+ ,4162.34
+ ,0
+ ,4136.22
+ ,0
+ ,4125.88
+ ,0
+ ,4031.48
+ ,0
+ ,3761.36
+ ,0
+ ,3408.56
+ ,0
+ ,3228.47
+ ,0
+ ,3090.45
+ ,0
+ ,2741.14
+ ,0
+ ,2980.44
+ ,0
+ ,3104.33
+ ,0
+ ,3181.57
+ ,0
+ ,2863.86
+ ,0
+ ,2898.01
+ ,0
+ ,3112.33
+ ,0
+ ,3254.33
+ ,0
+ ,3513.47
+ ,0
+ ,3587.61
+ ,0
+ ,3727.45
+ ,0
+ ,3793.34
+ ,0
+ ,3817.58
+ ,0
+ ,3845.13
+ ,0
+ ,3931.86
+ ,0
+ ,4197.52
+ ,0
+ ,4307.13
+ ,0
+ ,4229.43
+ ,0
+ ,4362.28
+ ,0
+ ,4217.34
+ ,0
+ ,4361.28
+ ,0
+ ,4327.74
+ ,0
+ ,4417.65
+ ,0
+ ,4557.68
+ ,0
+ ,4650.35
+ ,0
+ ,4967.18
+ ,0
+ ,5123.42
+ ,0
+ ,5290.85
+ ,0
+ ,5535.66
+ ,0
+ ,5514.06
+ ,0
+ ,5493.88
+ ,0
+ ,5694.83
+ ,0
+ ,5850.41
+ ,0
+ ,6116.64
+ ,0
+ ,6175.00
+ ,0
+ ,6513.58
+ ,0
+ ,6383.78
+ ,0
+ ,6673.66
+ ,0
+ ,6936.61
+ ,0
+ ,7300.68
+ ,0
+ ,7392.93
+ ,0
+ ,7497.31
+ ,0
+ ,7584.71
+ ,0
+ ,7160.79
+ ,0
+ ,7196.19
+ ,0
+ ,7245.63
+ ,0
+ ,7347.51
+ ,0
+ ,7425.75
+ ,0
+ ,7778.51
+ ,0
+ ,7822.33
+ ,0
+ ,8181.22
+ ,0
+ ,8371.47
+ ,0
+ ,8347.71
+ ,0
+ ,8672.11
+ ,0
+ ,8802.79
+ ,0
+ ,9138.46
+ ,0
+ ,9123.29
+ ,0
+ ,9023.21
+ ,1
+ ,8850.41
+ ,1
+ ,8864.58
+ ,1
+ ,9163.74
+ ,1
+ ,8516.66
+ ,1
+ ,8553.44
+ ,1
+ ,7555.20
+ ,1
+ ,7851.22
+ ,1
+ ,7442.00
+ ,1
+ ,7992.53
+ ,1
+ ,8264.04
+ ,1
+ ,7517.39
+ ,1
+ ,7200.40
+ ,1
+ ,7193.69
+ ,1
+ ,6193.58
+ ,1
+ ,5104.21
+ ,1
+ ,4800.46
+ ,1
+ ,4461.61
+ ,1
+ ,4398.59
+ ,1
+ ,4243.63
+ ,1
+ ,4293.82
+ ,1)
+ ,dim=c(2
+ ,108)
+ ,dimnames=list(c('Y'
+ ,'X')
+ ,1:108))
> y <- array(NA,dim=c(2,108),dimnames=list(c('Y','X'),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
Y X
1 4716.99 0
2 4926.65 0
3 4920.10 0
4 5170.09 0
5 5246.24 0
6 5283.61 0
7 4979.05 0
8 4825.20 0
9 4695.12 0
10 4711.54 0
11 4727.22 0
12 4384.96 0
13 4378.75 0
14 4472.93 0
15 4564.07 0
16 4310.54 0
17 4171.38 0
18 4049.38 0
19 3591.37 0
20 3720.46 0
21 4107.23 0
22 4101.71 0
23 4162.34 0
24 4136.22 0
25 4125.88 0
26 4031.48 0
27 3761.36 0
28 3408.56 0
29 3228.47 0
30 3090.45 0
31 2741.14 0
32 2980.44 0
33 3104.33 0
34 3181.57 0
35 2863.86 0
36 2898.01 0
37 3112.33 0
38 3254.33 0
39 3513.47 0
40 3587.61 0
41 3727.45 0
42 3793.34 0
43 3817.58 0
44 3845.13 0
45 3931.86 0
46 4197.52 0
47 4307.13 0
48 4229.43 0
49 4362.28 0
50 4217.34 0
51 4361.28 0
52 4327.74 0
53 4417.65 0
54 4557.68 0
55 4650.35 0
56 4967.18 0
57 5123.42 0
58 5290.85 0
59 5535.66 0
60 5514.06 0
61 5493.88 0
62 5694.83 0
63 5850.41 0
64 6116.64 0
65 6175.00 0
66 6513.58 0
67 6383.78 0
68 6673.66 0
69 6936.61 0
70 7300.68 0
71 7392.93 0
72 7497.31 0
73 7584.71 0
74 7160.79 0
75 7196.19 0
76 7245.63 0
77 7347.51 0
78 7425.75 0
79 7778.51 0
80 7822.33 0
81 8181.22 0
82 8371.47 0
83 8347.71 0
84 8672.11 0
85 8802.79 0
86 9138.46 0
87 9123.29 0
88 9023.21 1
89 8850.41 1
90 8864.58 1
91 9163.74 1
92 8516.66 1
93 8553.44 1
94 7555.20 1
95 7851.22 1
96 7442.00 1
97 7992.53 1
98 8264.04 1
99 7517.39 1
100 7200.40 1
101 7193.69 1
102 6193.58 1
103 5104.21 1
104 4800.46 1
105 4461.61 1
106 4398.59 1
107 4243.63 1
108 4293.82 1
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X
5157 1866
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-2779.4 -1150.2 -434.7 1269.9 3981.7
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5156.8 182.0 28.329 < 2e-16 ***
X 1866.3 412.8 4.521 1.61e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1698 on 106 degrees of freedom
Multiple R-squared: 0.1617, Adjusted R-squared: 0.1537
F-statistic: 20.44 on 1 and 106 DF, p-value: 1.609e-05
> 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.865875e-03 7.731751e-03 0.996134125
[2,] 8.331926e-04 1.666385e-03 0.999166807
[3,] 1.042054e-04 2.084108e-04 0.999895795
[4,] 1.764088e-05 3.528176e-05 0.999982359
[5,] 4.832997e-06 9.665994e-06 0.999995167
[6,] 1.037307e-06 2.074614e-06 0.999998963
[7,] 1.901765e-07 3.803531e-07 0.999999810
[8,] 2.004687e-07 4.009374e-07 0.999999800
[9,] 1.221503e-07 2.443006e-07 0.999999878
[10,] 3.868197e-08 7.736394e-08 0.999999961
[11,] 8.521867e-09 1.704373e-08 0.999999991
[12,] 4.419740e-09 8.839480e-09 0.999999996
[13,] 3.534118e-09 7.068237e-09 0.999999996
[14,] 3.691870e-09 7.383740e-09 0.999999996
[15,] 2.354866e-08 4.709731e-08 0.999999976
[16,] 3.518001e-08 7.036003e-08 0.999999965
[17,] 1.445025e-08 2.890051e-08 0.999999986
[18,] 5.765533e-09 1.153107e-08 0.999999994
[19,] 2.001165e-09 4.002329e-09 0.999999998
[20,] 7.059264e-10 1.411853e-09 0.999999999
[21,] 2.469504e-10 4.939008e-10 1.000000000
[22,] 1.004140e-10 2.008280e-10 1.000000000
[23,] 7.731142e-11 1.546228e-10 1.000000000
[24,] 1.714072e-10 3.428143e-10 1.000000000
[25,] 5.327347e-10 1.065469e-09 0.999999999
[26,] 1.854926e-09 3.709853e-09 0.999999998
[27,] 1.408707e-08 2.817413e-08 0.999999986
[28,] 3.352774e-08 6.705547e-08 0.999999966
[29,] 5.158675e-08 1.031735e-07 0.999999948
[30,] 6.392561e-08 1.278512e-07 0.999999936
[31,] 1.456270e-07 2.912540e-07 0.999999854
[32,] 2.840624e-07 5.681247e-07 0.999999716
[33,] 3.653611e-07 7.307223e-07 0.999999635
[34,] 3.856218e-07 7.712436e-07 0.999999614
[35,] 3.067626e-07 6.135252e-07 0.999999693
[36,] 2.353091e-07 4.706183e-07 0.999999765
[37,] 1.665123e-07 3.330247e-07 0.999999833
[38,] 1.169009e-07 2.338019e-07 0.999999883
[39,] 8.436512e-08 1.687302e-07 0.999999916
[40,] 6.282652e-08 1.256530e-07 0.999999937
[41,] 4.681158e-08 9.362316e-08 0.999999953
[42,] 3.325349e-08 6.650697e-08 0.999999967
[43,] 2.445250e-08 4.890501e-08 0.999999976
[44,] 1.908307e-08 3.816613e-08 0.999999981
[45,] 1.557599e-08 3.115197e-08 0.999999984
[46,] 1.392110e-08 2.784220e-08 0.999999986
[47,] 1.309116e-08 2.618232e-08 0.999999987
[48,] 1.351305e-08 2.702610e-08 0.999999986
[49,] 1.517589e-08 3.035177e-08 0.999999985
[50,] 1.877343e-08 3.754685e-08 0.999999981
[51,] 2.591136e-08 5.182271e-08 0.999999974
[52,] 4.388883e-08 8.777765e-08 0.999999956
[53,] 8.631133e-08 1.726227e-07 0.999999914
[54,] 1.971841e-07 3.943682e-07 0.999999803
[55,] 5.535262e-07 1.107052e-06 0.999999446
[56,] 1.421128e-06 2.842257e-06 0.999998579
[57,] 3.474754e-06 6.949509e-06 0.999996525
[58,] 9.446360e-06 1.889272e-05 0.999990554
[59,] 2.670728e-05 5.341455e-05 0.999973293
[60,] 8.234181e-05 1.646836e-04 0.999917658
[61,] 2.270437e-04 4.540873e-04 0.999772956
[62,] 6.727167e-04 1.345433e-03 0.999327283
[63,] 1.534613e-03 3.069226e-03 0.998465387
[64,] 3.526145e-03 7.052290e-03 0.996473855
[65,] 7.803058e-03 1.560612e-02 0.992196942
[66,] 1.713643e-02 3.427285e-02 0.982863574
[67,] 3.185808e-02 6.371615e-02 0.968141923
[68,] 5.230359e-02 1.046072e-01 0.947696406
[69,] 7.742082e-02 1.548416e-01 0.922579178
[70,] 9.560971e-02 1.912194e-01 0.904390291
[71,] 1.143199e-01 2.286399e-01 0.885680067
[72,] 1.334747e-01 2.669494e-01 0.866525281
[73,] 1.533916e-01 3.067832e-01 0.846608412
[74,] 1.736326e-01 3.472652e-01 0.826367419
[75,] 1.969627e-01 3.939253e-01 0.803037348
[76,] 2.181432e-01 4.362864e-01 0.781856806
[77,] 2.423722e-01 4.847444e-01 0.757627778
[78,] 2.654280e-01 5.308559e-01 0.734572032
[79,] 2.815938e-01 5.631875e-01 0.718406236
[80,] 2.994291e-01 5.988583e-01 0.700570874
[81,] 3.131758e-01 6.263516e-01 0.686824218
[82,] 3.291195e-01 6.582390e-01 0.670880514
[83,] 3.348593e-01 6.697187e-01 0.665140667
[84,] 3.408177e-01 6.816354e-01 0.659182293
[85,] 3.436587e-01 6.873175e-01 0.656341253
[86,] 3.582291e-01 7.164582e-01 0.641770888
[87,] 4.211015e-01 8.422030e-01 0.578898515
[88,] 4.393923e-01 8.787845e-01 0.560607739
[89,] 4.809239e-01 9.618478e-01 0.519076097
[90,] 4.462855e-01 8.925710e-01 0.553714500
[91,] 4.422065e-01 8.844130e-01 0.557793494
[92,] 4.131912e-01 8.263824e-01 0.586808824
[93,] 4.601276e-01 9.202552e-01 0.539872394
[94,] 6.134265e-01 7.731471e-01 0.386573548
[95,] 7.000716e-01 5.998567e-01 0.299928370
[96,] 7.909802e-01 4.180396e-01 0.209019789
[97,] 9.497352e-01 1.005296e-01 0.050264789
[98,] 9.942385e-01 1.152293e-02 0.005761465
[99,] 9.947121e-01 1.057579e-02 0.005287893
> postscript(file="/var/www/html/rcomp/tmp/1aqa71258644287.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/2x3101258644287.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/34qm01258644287.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/46n3c1258644287.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/5t6vq1258644287.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
-439.78184 -230.12184 -236.67184 13.31816 89.46816 126.83816
7 8 9 10 11 12
-177.72184 -331.57184 -461.65184 -445.23184 -429.55184 -771.81184
13 14 15 16 17 18
-778.02184 -683.84184 -592.70184 -846.23184 -985.39184 -1107.39184
19 20 21 22 23 24
-1565.40184 -1436.31184 -1049.54184 -1055.06184 -994.43184 -1020.55184
25 26 27 28 29 30
-1030.89184 -1125.29184 -1395.41184 -1748.21184 -1928.30184 -2066.32184
31 32 33 34 35 36
-2415.63184 -2176.33184 -2052.44184 -1975.20184 -2292.91184 -2258.76184
37 38 39 40 41 42
-2044.44184 -1902.44184 -1643.30184 -1569.16184 -1429.32184 -1363.43184
43 44 45 46 47 48
-1339.19184 -1311.64184 -1224.91184 -959.25184 -849.64184 -927.34184
49 50 51 52 53 54
-794.49184 -939.43184 -795.49184 -829.03184 -739.12184 -599.09184
55 56 57 58 59 60
-506.42184 -189.59184 -33.35184 134.07816 378.88816 357.28816
61 62 63 64 65 66
337.10816 538.05816 693.63816 959.86816 1018.22816 1356.80816
67 68 69 70 71 72
1227.00816 1516.88816 1779.83816 2143.90816 2236.15816 2340.53816
73 74 75 76 77 78
2427.93816 2004.01816 2039.41816 2088.85816 2190.73816 2268.97816
79 80 81 82 83 84
2621.73816 2665.55816 3024.44816 3214.69816 3190.93816 3515.33816
85 86 87 88 89 90
3646.01816 3981.68816 3966.51816 2000.14286 1827.34286 1841.51286
91 92 93 94 95 96
2140.67286 1493.59286 1530.37286 532.13286 828.15286 418.93286
97 98 99 100 101 102
969.46286 1240.97286 494.32286 177.33286 170.62286 -829.48714
103 104 105 106 107 108
-1918.85714 -2222.60714 -2561.45714 -2624.47714 -2779.43714 -2729.24714
> postscript(file="/var/www/html/rcomp/tmp/63h5c1258644287.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 -439.78184 NA
1 -230.12184 -439.78184
2 -236.67184 -230.12184
3 13.31816 -236.67184
4 89.46816 13.31816
5 126.83816 89.46816
6 -177.72184 126.83816
7 -331.57184 -177.72184
8 -461.65184 -331.57184
9 -445.23184 -461.65184
10 -429.55184 -445.23184
11 -771.81184 -429.55184
12 -778.02184 -771.81184
13 -683.84184 -778.02184
14 -592.70184 -683.84184
15 -846.23184 -592.70184
16 -985.39184 -846.23184
17 -1107.39184 -985.39184
18 -1565.40184 -1107.39184
19 -1436.31184 -1565.40184
20 -1049.54184 -1436.31184
21 -1055.06184 -1049.54184
22 -994.43184 -1055.06184
23 -1020.55184 -994.43184
24 -1030.89184 -1020.55184
25 -1125.29184 -1030.89184
26 -1395.41184 -1125.29184
27 -1748.21184 -1395.41184
28 -1928.30184 -1748.21184
29 -2066.32184 -1928.30184
30 -2415.63184 -2066.32184
31 -2176.33184 -2415.63184
32 -2052.44184 -2176.33184
33 -1975.20184 -2052.44184
34 -2292.91184 -1975.20184
35 -2258.76184 -2292.91184
36 -2044.44184 -2258.76184
37 -1902.44184 -2044.44184
38 -1643.30184 -1902.44184
39 -1569.16184 -1643.30184
40 -1429.32184 -1569.16184
41 -1363.43184 -1429.32184
42 -1339.19184 -1363.43184
43 -1311.64184 -1339.19184
44 -1224.91184 -1311.64184
45 -959.25184 -1224.91184
46 -849.64184 -959.25184
47 -927.34184 -849.64184
48 -794.49184 -927.34184
49 -939.43184 -794.49184
50 -795.49184 -939.43184
51 -829.03184 -795.49184
52 -739.12184 -829.03184
53 -599.09184 -739.12184
54 -506.42184 -599.09184
55 -189.59184 -506.42184
56 -33.35184 -189.59184
57 134.07816 -33.35184
58 378.88816 134.07816
59 357.28816 378.88816
60 337.10816 357.28816
61 538.05816 337.10816
62 693.63816 538.05816
63 959.86816 693.63816
64 1018.22816 959.86816
65 1356.80816 1018.22816
66 1227.00816 1356.80816
67 1516.88816 1227.00816
68 1779.83816 1516.88816
69 2143.90816 1779.83816
70 2236.15816 2143.90816
71 2340.53816 2236.15816
72 2427.93816 2340.53816
73 2004.01816 2427.93816
74 2039.41816 2004.01816
75 2088.85816 2039.41816
76 2190.73816 2088.85816
77 2268.97816 2190.73816
78 2621.73816 2268.97816
79 2665.55816 2621.73816
80 3024.44816 2665.55816
81 3214.69816 3024.44816
82 3190.93816 3214.69816
83 3515.33816 3190.93816
84 3646.01816 3515.33816
85 3981.68816 3646.01816
86 3966.51816 3981.68816
87 2000.14286 3966.51816
88 1827.34286 2000.14286
89 1841.51286 1827.34286
90 2140.67286 1841.51286
91 1493.59286 2140.67286
92 1530.37286 1493.59286
93 532.13286 1530.37286
94 828.15286 532.13286
95 418.93286 828.15286
96 969.46286 418.93286
97 1240.97286 969.46286
98 494.32286 1240.97286
99 177.33286 494.32286
100 170.62286 177.33286
101 -829.48714 170.62286
102 -1918.85714 -829.48714
103 -2222.60714 -1918.85714
104 -2561.45714 -2222.60714
105 -2624.47714 -2561.45714
106 -2779.43714 -2624.47714
107 -2729.24714 -2779.43714
108 NA -2729.24714
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -230.12184 -439.78184
[2,] -236.67184 -230.12184
[3,] 13.31816 -236.67184
[4,] 89.46816 13.31816
[5,] 126.83816 89.46816
[6,] -177.72184 126.83816
[7,] -331.57184 -177.72184
[8,] -461.65184 -331.57184
[9,] -445.23184 -461.65184
[10,] -429.55184 -445.23184
[11,] -771.81184 -429.55184
[12,] -778.02184 -771.81184
[13,] -683.84184 -778.02184
[14,] -592.70184 -683.84184
[15,] -846.23184 -592.70184
[16,] -985.39184 -846.23184
[17,] -1107.39184 -985.39184
[18,] -1565.40184 -1107.39184
[19,] -1436.31184 -1565.40184
[20,] -1049.54184 -1436.31184
[21,] -1055.06184 -1049.54184
[22,] -994.43184 -1055.06184
[23,] -1020.55184 -994.43184
[24,] -1030.89184 -1020.55184
[25,] -1125.29184 -1030.89184
[26,] -1395.41184 -1125.29184
[27,] -1748.21184 -1395.41184
[28,] -1928.30184 -1748.21184
[29,] -2066.32184 -1928.30184
[30,] -2415.63184 -2066.32184
[31,] -2176.33184 -2415.63184
[32,] -2052.44184 -2176.33184
[33,] -1975.20184 -2052.44184
[34,] -2292.91184 -1975.20184
[35,] -2258.76184 -2292.91184
[36,] -2044.44184 -2258.76184
[37,] -1902.44184 -2044.44184
[38,] -1643.30184 -1902.44184
[39,] -1569.16184 -1643.30184
[40,] -1429.32184 -1569.16184
[41,] -1363.43184 -1429.32184
[42,] -1339.19184 -1363.43184
[43,] -1311.64184 -1339.19184
[44,] -1224.91184 -1311.64184
[45,] -959.25184 -1224.91184
[46,] -849.64184 -959.25184
[47,] -927.34184 -849.64184
[48,] -794.49184 -927.34184
[49,] -939.43184 -794.49184
[50,] -795.49184 -939.43184
[51,] -829.03184 -795.49184
[52,] -739.12184 -829.03184
[53,] -599.09184 -739.12184
[54,] -506.42184 -599.09184
[55,] -189.59184 -506.42184
[56,] -33.35184 -189.59184
[57,] 134.07816 -33.35184
[58,] 378.88816 134.07816
[59,] 357.28816 378.88816
[60,] 337.10816 357.28816
[61,] 538.05816 337.10816
[62,] 693.63816 538.05816
[63,] 959.86816 693.63816
[64,] 1018.22816 959.86816
[65,] 1356.80816 1018.22816
[66,] 1227.00816 1356.80816
[67,] 1516.88816 1227.00816
[68,] 1779.83816 1516.88816
[69,] 2143.90816 1779.83816
[70,] 2236.15816 2143.90816
[71,] 2340.53816 2236.15816
[72,] 2427.93816 2340.53816
[73,] 2004.01816 2427.93816
[74,] 2039.41816 2004.01816
[75,] 2088.85816 2039.41816
[76,] 2190.73816 2088.85816
[77,] 2268.97816 2190.73816
[78,] 2621.73816 2268.97816
[79,] 2665.55816 2621.73816
[80,] 3024.44816 2665.55816
[81,] 3214.69816 3024.44816
[82,] 3190.93816 3214.69816
[83,] 3515.33816 3190.93816
[84,] 3646.01816 3515.33816
[85,] 3981.68816 3646.01816
[86,] 3966.51816 3981.68816
[87,] 2000.14286 3966.51816
[88,] 1827.34286 2000.14286
[89,] 1841.51286 1827.34286
[90,] 2140.67286 1841.51286
[91,] 1493.59286 2140.67286
[92,] 1530.37286 1493.59286
[93,] 532.13286 1530.37286
[94,] 828.15286 532.13286
[95,] 418.93286 828.15286
[96,] 969.46286 418.93286
[97,] 1240.97286 969.46286
[98,] 494.32286 1240.97286
[99,] 177.33286 494.32286
[100,] 170.62286 177.33286
[101,] -829.48714 170.62286
[102,] -1918.85714 -829.48714
[103,] -2222.60714 -1918.85714
[104,] -2561.45714 -2222.60714
[105,] -2624.47714 -2561.45714
[106,] -2779.43714 -2624.47714
[107,] -2729.24714 -2779.43714
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -230.12184 -439.78184
2 -236.67184 -230.12184
3 13.31816 -236.67184
4 89.46816 13.31816
5 126.83816 89.46816
6 -177.72184 126.83816
7 -331.57184 -177.72184
8 -461.65184 -331.57184
9 -445.23184 -461.65184
10 -429.55184 -445.23184
11 -771.81184 -429.55184
12 -778.02184 -771.81184
13 -683.84184 -778.02184
14 -592.70184 -683.84184
15 -846.23184 -592.70184
16 -985.39184 -846.23184
17 -1107.39184 -985.39184
18 -1565.40184 -1107.39184
19 -1436.31184 -1565.40184
20 -1049.54184 -1436.31184
21 -1055.06184 -1049.54184
22 -994.43184 -1055.06184
23 -1020.55184 -994.43184
24 -1030.89184 -1020.55184
25 -1125.29184 -1030.89184
26 -1395.41184 -1125.29184
27 -1748.21184 -1395.41184
28 -1928.30184 -1748.21184
29 -2066.32184 -1928.30184
30 -2415.63184 -2066.32184
31 -2176.33184 -2415.63184
32 -2052.44184 -2176.33184
33 -1975.20184 -2052.44184
34 -2292.91184 -1975.20184
35 -2258.76184 -2292.91184
36 -2044.44184 -2258.76184
37 -1902.44184 -2044.44184
38 -1643.30184 -1902.44184
39 -1569.16184 -1643.30184
40 -1429.32184 -1569.16184
41 -1363.43184 -1429.32184
42 -1339.19184 -1363.43184
43 -1311.64184 -1339.19184
44 -1224.91184 -1311.64184
45 -959.25184 -1224.91184
46 -849.64184 -959.25184
47 -927.34184 -849.64184
48 -794.49184 -927.34184
49 -939.43184 -794.49184
50 -795.49184 -939.43184
51 -829.03184 -795.49184
52 -739.12184 -829.03184
53 -599.09184 -739.12184
54 -506.42184 -599.09184
55 -189.59184 -506.42184
56 -33.35184 -189.59184
57 134.07816 -33.35184
58 378.88816 134.07816
59 357.28816 378.88816
60 337.10816 357.28816
61 538.05816 337.10816
62 693.63816 538.05816
63 959.86816 693.63816
64 1018.22816 959.86816
65 1356.80816 1018.22816
66 1227.00816 1356.80816
67 1516.88816 1227.00816
68 1779.83816 1516.88816
69 2143.90816 1779.83816
70 2236.15816 2143.90816
71 2340.53816 2236.15816
72 2427.93816 2340.53816
73 2004.01816 2427.93816
74 2039.41816 2004.01816
75 2088.85816 2039.41816
76 2190.73816 2088.85816
77 2268.97816 2190.73816
78 2621.73816 2268.97816
79 2665.55816 2621.73816
80 3024.44816 2665.55816
81 3214.69816 3024.44816
82 3190.93816 3214.69816
83 3515.33816 3190.93816
84 3646.01816 3515.33816
85 3981.68816 3646.01816
86 3966.51816 3981.68816
87 2000.14286 3966.51816
88 1827.34286 2000.14286
89 1841.51286 1827.34286
90 2140.67286 1841.51286
91 1493.59286 2140.67286
92 1530.37286 1493.59286
93 532.13286 1530.37286
94 828.15286 532.13286
95 418.93286 828.15286
96 969.46286 418.93286
97 1240.97286 969.46286
98 494.32286 1240.97286
99 177.33286 494.32286
100 170.62286 177.33286
101 -829.48714 170.62286
102 -1918.85714 -829.48714
103 -2222.60714 -1918.85714
104 -2561.45714 -2222.60714
105 -2624.47714 -2561.45714
106 -2779.43714 -2624.47714
107 -2729.24714 -2779.43714
> 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/7vj7b1258644287.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/89k661258644287.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/9w2si1258644287.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/10qbbw1258644287.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/11l31x1258644287.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/12vsax1258644287.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/13y1kr1258644287.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/14w5xi1258644287.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/15yssj1258644287.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/16fcdr1258644287.tab")
+ }
>
> system("convert tmp/1aqa71258644287.ps tmp/1aqa71258644287.png")
> system("convert tmp/2x3101258644287.ps tmp/2x3101258644287.png")
> system("convert tmp/34qm01258644287.ps tmp/34qm01258644287.png")
> system("convert tmp/46n3c1258644287.ps tmp/46n3c1258644287.png")
> system("convert tmp/5t6vq1258644287.ps tmp/5t6vq1258644287.png")
> system("convert tmp/63h5c1258644287.ps tmp/63h5c1258644287.png")
> system("convert tmp/7vj7b1258644287.ps tmp/7vj7b1258644287.png")
> system("convert tmp/89k661258644287.ps tmp/89k661258644287.png")
> system("convert tmp/9w2si1258644287.ps tmp/9w2si1258644287.png")
> system("convert tmp/10qbbw1258644287.ps tmp/10qbbw1258644287.png")
>
>
> proc.time()
user system elapsed
3.050 1.635 4.051