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(106370
+ ,100.3
+ ,109375
+ ,101.9
+ ,116476
+ ,102.1
+ ,123297
+ ,103.2
+ ,114813
+ ,103.7
+ ,117925
+ ,106.2
+ ,126466
+ ,107.7
+ ,131235
+ ,109.9
+ ,120546
+ ,111.7
+ ,123791
+ ,114.9
+ ,129813
+ ,116
+ ,133463
+ ,118.3
+ ,122987
+ ,120.4
+ ,125418
+ ,126
+ ,130199
+ ,128.1
+ ,133016
+ ,130.1
+ ,121454
+ ,130.8
+ ,122044
+ ,133.6
+ ,128313
+ ,134.2
+ ,131556
+ ,135.5
+ ,120027
+ ,136.2
+ ,123001
+ ,139.1
+ ,130111
+ ,139
+ ,132524
+ ,139.6
+ ,123742
+ ,138.7
+ ,124931
+ ,140.9
+ ,133646
+ ,141.3
+ ,136557
+ ,141.8
+ ,127509
+ ,142
+ ,128945
+ ,144.5
+ ,137191
+ ,144.6
+ ,139716
+ ,145.5
+ ,129083
+ ,146.8
+ ,131604
+ ,149.5
+ ,139413
+ ,149.9
+ ,143125
+ ,150.1
+ ,133948
+ ,150.9
+ ,137116
+ ,152.8
+ ,144864
+ ,153.1
+ ,149277
+ ,154
+ ,138796
+ ,154.9
+ ,143258
+ ,156.9
+ ,150034
+ ,158.4
+ ,154708
+ ,159.7
+ ,144888
+ ,160.2
+ ,148762
+ ,163.2
+ ,156500
+ ,163.7
+ ,161088
+ ,164.4
+ ,152772
+ ,163.7
+ ,158011
+ ,165.5
+ ,163318
+ ,165.6
+ ,169969
+ ,166.8
+ ,162269
+ ,167.5
+ ,165765
+ ,170.6
+ ,170600
+ ,170.9
+ ,174681
+ ,172
+ ,166364
+ ,171.8
+ ,170240
+ ,173.9
+ ,176150
+ ,174
+ ,182056
+ ,173.8
+ ,172218
+ ,173.9
+ ,177856
+ ,176
+ ,182253
+ ,176.6
+ ,188090
+ ,178.2
+ ,176863
+ ,179.2
+ ,183273
+ ,181.3
+ ,187969
+ ,181.8
+ ,194650
+ ,182.9
+ ,183036
+ ,183.8
+ ,189516
+ ,186.3
+ ,193805
+ ,187.4
+ ,200499
+ ,189.2
+ ,188142
+ ,189.7
+ ,193732
+ ,191.9
+ ,197126
+ ,192.6
+ ,205140
+ ,193.7
+ ,191751
+ ,194.2
+ ,196700
+ ,197.6
+ ,199784
+ ,199.3
+ ,207360
+ ,201.4
+ ,196101
+ ,203
+ ,200824
+ ,206.3
+ ,205743
+ ,207.1
+ ,212489
+ ,209.8
+ ,200810
+ ,211.1
+ ,203683
+ ,215.3
+ ,207286
+ ,217.4
+ ,210910
+ ,215.5
+ ,194915
+ ,210.9
+ ,217920
+ ,212.6)
+ ,dim=c(2
+ ,90)
+ ,dimnames=list(c('HFCE'
+ ,'RPI')
+ ,1:90))
> y <- array(NA,dim=c(2,90),dimnames=list(c('HFCE','RPI'),1:90))
> 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 = 'Include Quarterly 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
HFCE RPI Q1 Q2 Q3
1 106370 100.3 1 0 0
2 109375 101.9 0 1 0
3 116476 102.1 0 0 1
4 123297 103.2 0 0 0
5 114813 103.7 1 0 0
6 117925 106.2 0 1 0
7 126466 107.7 0 0 1
8 131235 109.9 0 0 0
9 120546 111.7 1 0 0
10 123791 114.9 0 1 0
11 129813 116.0 0 0 1
12 133463 118.3 0 0 0
13 122987 120.4 1 0 0
14 125418 126.0 0 1 0
15 130199 128.1 0 0 1
16 133016 130.1 0 0 0
17 121454 130.8 1 0 0
18 122044 133.6 0 1 0
19 128313 134.2 0 0 1
20 131556 135.5 0 0 0
21 120027 136.2 1 0 0
22 123001 139.1 0 1 0
23 130111 139.0 0 0 1
24 132524 139.6 0 0 0
25 123742 138.7 1 0 0
26 124931 140.9 0 1 0
27 133646 141.3 0 0 1
28 136557 141.8 0 0 0
29 127509 142.0 1 0 0
30 128945 144.5 0 1 0
31 137191 144.6 0 0 1
32 139716 145.5 0 0 0
33 129083 146.8 1 0 0
34 131604 149.5 0 1 0
35 139413 149.9 0 0 1
36 143125 150.1 0 0 0
37 133948 150.9 1 0 0
38 137116 152.8 0 1 0
39 144864 153.1 0 0 1
40 149277 154.0 0 0 0
41 138796 154.9 1 0 0
42 143258 156.9 0 1 0
43 150034 158.4 0 0 1
44 154708 159.7 0 0 0
45 144888 160.2 1 0 0
46 148762 163.2 0 1 0
47 156500 163.7 0 0 1
48 161088 164.4 0 0 0
49 152772 163.7 1 0 0
50 158011 165.5 0 1 0
51 163318 165.6 0 0 1
52 169969 166.8 0 0 0
53 162269 167.5 1 0 0
54 165765 170.6 0 1 0
55 170600 170.9 0 0 1
56 174681 172.0 0 0 0
57 166364 171.8 1 0 0
58 170240 173.9 0 1 0
59 176150 174.0 0 0 1
60 182056 173.8 0 0 0
61 172218 173.9 1 0 0
62 177856 176.0 0 1 0
63 182253 176.6 0 0 1
64 188090 178.2 0 0 0
65 176863 179.2 1 0 0
66 183273 181.3 0 1 0
67 187969 181.8 0 0 1
68 194650 182.9 0 0 0
69 183036 183.8 1 0 0
70 189516 186.3 0 1 0
71 193805 187.4 0 0 1
72 200499 189.2 0 0 0
73 188142 189.7 1 0 0
74 193732 191.9 0 1 0
75 197126 192.6 0 0 1
76 205140 193.7 0 0 0
77 191751 194.2 1 0 0
78 196700 197.6 0 1 0
79 199784 199.3 0 0 1
80 207360 201.4 0 0 0
81 196101 203.0 1 0 0
82 200824 206.3 0 1 0
83 205743 207.1 0 0 1
84 212489 209.8 0 0 0
85 200810 211.1 1 0 0
86 203683 215.3 0 1 0
87 207286 217.4 0 0 1
88 210910 215.5 0 0 0
89 194915 210.9 1 0 0
90 217920 212.6 0 1 0
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) RPI Q1 Q2 Q3
13234.1 940.3 -10550.7 -8526.7 -3873.4
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-13674.9 -8545.2 596.1 7471.9 15837.5
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 13234.08 5055.52 2.618 0.010477 *
RPI 940.28 29.26 32.130 < 2e-16 ***
Q1 -10550.72 2595.90 -4.064 0.000107 ***
Q2 -8526.68 2595.13 -3.286 0.001480 **
Q3 -3873.42 2623.98 -1.476 0.143595
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 8702 on 85 degrees of freedom
Multiple R-squared: 0.9256, Adjusted R-squared: 0.9221
F-statistic: 264.4 on 4 and 85 DF, p-value: < 2.2e-16
> if (n > n25) {
+ kp3 <- k + 3
+ nmkm3 <- n - k - 3
+ gqarr <- array(NA, dim=c(nmkm3-kp3+1,3))
+ numgqtests <- 0
+ numsignificant1 <- 0
+ numsignificant5 <- 0
+ numsignificant10 <- 0
+ for (mypoint in kp3:nmkm3) {
+ j <- 0
+ numgqtests <- numgqtests + 1
+ for (myalt in c('greater', 'two.sided', 'less')) {
+ j <- j + 1
+ gqarr[mypoint-kp3+1,j] <- gqtest(mylm, point=mypoint, alternative=myalt)$p.value
+ }
+ if (gqarr[mypoint-kp3+1,2] < 0.01) numsignificant1 <- numsignificant1 + 1
+ if (gqarr[mypoint-kp3+1,2] < 0.05) numsignificant5 <- numsignificant5 + 1
+ if (gqarr[mypoint-kp3+1,2] < 0.10) numsignificant10 <- numsignificant10 + 1
+ }
+ gqarr
+ }
[,1] [,2] [,3]
[1,] 0.0192848176 3.856964e-02 9.807152e-01
[2,] 0.0181763932 3.635279e-02 9.818236e-01
[3,] 0.0112702234 2.254045e-02 9.887298e-01
[4,] 0.0102077111 2.041542e-02 9.897923e-01
[5,] 0.0158809619 3.176192e-02 9.841190e-01
[6,] 0.0172090380 3.441808e-02 9.827910e-01
[7,] 0.0251250166 5.025003e-02 9.748750e-01
[8,] 0.0455787233 9.115745e-02 9.544213e-01
[9,] 0.0735902740 1.471805e-01 9.264097e-01
[10,] 0.0751607088 1.503214e-01 9.248393e-01
[11,] 0.0834701145 1.669402e-01 9.165299e-01
[12,] 0.0798232834 1.596466e-01 9.201767e-01
[13,] 0.0778766171 1.557532e-01 9.221234e-01
[14,] 0.0650559044 1.301118e-01 9.349441e-01
[15,] 0.0496304190 9.926084e-02 9.503696e-01
[16,] 0.0335536142 6.710723e-02 9.664464e-01
[17,] 0.0253091709 5.061834e-02 9.746908e-01
[18,] 0.0155650767 3.113015e-02 9.844349e-01
[19,] 0.0097555519 1.951110e-02 9.902444e-01
[20,] 0.0059315635 1.186313e-02 9.940684e-01
[21,] 0.0034457179 6.891436e-03 9.965543e-01
[22,] 0.0023530490 4.706098e-03 9.976470e-01
[23,] 0.0015384688 3.076938e-03 9.984615e-01
[24,] 0.0010948560 2.189712e-03 9.989051e-01
[25,] 0.0006920751 1.384150e-03 9.993079e-01
[26,] 0.0004597438 9.194875e-04 9.995403e-01
[27,] 0.0003983694 7.967387e-04 9.996016e-01
[28,] 0.0003091724 6.183449e-04 9.996908e-01
[29,] 0.0002755501 5.511001e-04 9.997244e-01
[30,] 0.0003209038 6.418077e-04 9.996791e-01
[31,] 0.0005706025 1.141205e-03 9.994294e-01
[32,] 0.0008237697 1.647539e-03 9.991762e-01
[33,] 0.0014508866 2.901773e-03 9.985491e-01
[34,] 0.0024468323 4.893665e-03 9.975532e-01
[35,] 0.0069514448 1.390289e-02 9.930486e-01
[36,] 0.0120405041 2.408101e-02 9.879595e-01
[37,] 0.0256840245 5.136805e-02 9.743160e-01
[38,] 0.0516792911 1.033586e-01 9.483207e-01
[39,] 0.1470940420 2.941881e-01 8.529060e-01
[40,] 0.2509792135 5.019584e-01 7.490208e-01
[41,] 0.4418843466 8.837687e-01 5.581157e-01
[42,] 0.6192918862 7.614162e-01 3.807081e-01
[43,] 0.8259142594 3.481715e-01 1.740857e-01
[44,] 0.8981427172 2.037146e-01 1.018573e-01
[45,] 0.9535171655 9.296567e-02 4.648283e-02
[46,] 0.9756600946 4.867981e-02 2.433991e-02
[47,] 0.9920432795 1.591344e-02 7.956721e-03
[48,] 0.9955139860 8.972028e-03 4.486014e-03
[49,] 0.9986206732 2.758654e-03 1.379327e-03
[50,] 0.9991787890 1.642422e-03 8.212110e-04
[51,] 0.9998358596 3.282809e-04 1.641404e-04
[52,] 0.9998998621 2.002758e-04 1.001379e-04
[53,] 0.9999470923 1.058154e-04 5.290772e-05
[54,] 0.9999485516 1.028969e-04 5.144843e-05
[55,] 0.9999710412 5.791751e-05 2.895876e-05
[56,] 0.9999651408 6.971843e-05 3.485921e-05
[57,] 0.9999670818 6.583637e-05 3.291818e-05
[58,] 0.9999583742 8.325165e-05 4.162583e-05
[59,] 0.9999681610 6.367806e-05 3.183903e-05
[60,] 0.9999451494 1.097012e-04 5.485060e-05
[61,] 0.9999108238 1.783524e-04 8.917620e-05
[62,] 0.9998260500 3.479000e-04 1.739500e-04
[63,] 0.9997512269 4.975462e-04 2.487731e-04
[64,] 0.9994620984 1.075803e-03 5.379016e-04
[65,] 0.9988865848 2.226830e-03 1.113415e-03
[66,] 0.9975671707 4.865659e-03 2.432829e-03
[67,] 0.9959208423 8.158315e-03 4.079158e-03
[68,] 0.9912108237 1.757835e-02 8.789176e-03
[69,] 0.9831555802 3.368884e-02 1.684442e-02
[70,] 0.9665722157 6.685557e-02 3.342778e-02
[71,] 0.9480754559 1.038491e-01 5.192454e-02
[72,] 0.9029747794 1.940504e-01 9.702522e-02
[73,] 0.8244187574 3.511625e-01 1.755812e-01
[74,] 0.7019244536 5.961511e-01 2.980755e-01
[75,] 0.6873027413 6.253945e-01 3.126973e-01
> postscript(file="/var/www/html/rcomp/tmp/191ti1258724451.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/21lkh1258724451.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/328qj1258724451.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/409dh1258724451.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/5n8bg1258724451.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 = 90
Frequency = 1
1 2 3 4 5 6
9376.8300 8853.3439 11113.0273 13026.3051 14622.8873 13360.1517
7 8 9 10 11 12
15837.4748 14664.4476 12833.6694 11045.7397 11380.1736 8994.1187
13 14 15 16 17 18
7094.2574 2235.6623 388.8190 -2548.1527 -4217.6259 -8284.4447
19 20 21 22 23 24
-7232.8722 -9085.6499 -10722.1230 -12498.9695 -9948.2029 -11972.7865
25 26 27 28 29 30
-9357.8162 -12261.4686 -8575.8406 -10008.3965 -8693.7311 -11632.4667
31 32 33 34 35 36
-8133.7555 -10328.4223 -11633.0618 -13674.8529 -10895.2249 -11244.6976
37 38 39 40 41 42
-10623.1985 -11265.7678 -8453.1121 -8759.7788 -9536.3075 -8978.9045
43 44 45 46 47 48
-8266.5814 -8688.3591 -8427.7769 -9398.6511 -6784.0508 -6727.6622
49 50 51 52 53 54
-3834.7472 -2312.2887 -1752.5776 -103.3275 2089.1993 646.2973
55 56 57 58 59 60
545.9530 -280.7692 2141.0071 2018.3824 3181.0936 5401.7318
61 62 63 64 65 66
6020.4249 7659.8002 6839.3727 7298.5119 5681.9555 8093.3308
67 68 69 70 71 72
7665.9311 9439.2089 7529.6802 9634.9446 8236.3785 9364.4622
73 74 75 76 77 78
7088.0445 8585.3921 6667.9369 9774.2146 6465.7969 6193.8118
79 80 81 82 83 84
3026.0793 4754.0799 2541.3572 2137.3998 1650.9169 1984.7510
85 86 87 88 89 90
-365.8885 -3466.0954 -6490.9387 -4953.8293 -6072.8331 13309.6531
> postscript(file="/var/www/html/rcomp/tmp/6buxr1258724451.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 = 90
Frequency = 1
lag(myerror, k = 1) myerror
0 9376.8300 NA
1 8853.3439 9376.8300
2 11113.0273 8853.3439
3 13026.3051 11113.0273
4 14622.8873 13026.3051
5 13360.1517 14622.8873
6 15837.4748 13360.1517
7 14664.4476 15837.4748
8 12833.6694 14664.4476
9 11045.7397 12833.6694
10 11380.1736 11045.7397
11 8994.1187 11380.1736
12 7094.2574 8994.1187
13 2235.6623 7094.2574
14 388.8190 2235.6623
15 -2548.1527 388.8190
16 -4217.6259 -2548.1527
17 -8284.4447 -4217.6259
18 -7232.8722 -8284.4447
19 -9085.6499 -7232.8722
20 -10722.1230 -9085.6499
21 -12498.9695 -10722.1230
22 -9948.2029 -12498.9695
23 -11972.7865 -9948.2029
24 -9357.8162 -11972.7865
25 -12261.4686 -9357.8162
26 -8575.8406 -12261.4686
27 -10008.3965 -8575.8406
28 -8693.7311 -10008.3965
29 -11632.4667 -8693.7311
30 -8133.7555 -11632.4667
31 -10328.4223 -8133.7555
32 -11633.0618 -10328.4223
33 -13674.8529 -11633.0618
34 -10895.2249 -13674.8529
35 -11244.6976 -10895.2249
36 -10623.1985 -11244.6976
37 -11265.7678 -10623.1985
38 -8453.1121 -11265.7678
39 -8759.7788 -8453.1121
40 -9536.3075 -8759.7788
41 -8978.9045 -9536.3075
42 -8266.5814 -8978.9045
43 -8688.3591 -8266.5814
44 -8427.7769 -8688.3591
45 -9398.6511 -8427.7769
46 -6784.0508 -9398.6511
47 -6727.6622 -6784.0508
48 -3834.7472 -6727.6622
49 -2312.2887 -3834.7472
50 -1752.5776 -2312.2887
51 -103.3275 -1752.5776
52 2089.1993 -103.3275
53 646.2973 2089.1993
54 545.9530 646.2973
55 -280.7692 545.9530
56 2141.0071 -280.7692
57 2018.3824 2141.0071
58 3181.0936 2018.3824
59 5401.7318 3181.0936
60 6020.4249 5401.7318
61 7659.8002 6020.4249
62 6839.3727 7659.8002
63 7298.5119 6839.3727
64 5681.9555 7298.5119
65 8093.3308 5681.9555
66 7665.9311 8093.3308
67 9439.2089 7665.9311
68 7529.6802 9439.2089
69 9634.9446 7529.6802
70 8236.3785 9634.9446
71 9364.4622 8236.3785
72 7088.0445 9364.4622
73 8585.3921 7088.0445
74 6667.9369 8585.3921
75 9774.2146 6667.9369
76 6465.7969 9774.2146
77 6193.8118 6465.7969
78 3026.0793 6193.8118
79 4754.0799 3026.0793
80 2541.3572 4754.0799
81 2137.3998 2541.3572
82 1650.9169 2137.3998
83 1984.7510 1650.9169
84 -365.8885 1984.7510
85 -3466.0954 -365.8885
86 -6490.9387 -3466.0954
87 -4953.8293 -6490.9387
88 -6072.8331 -4953.8293
89 13309.6531 -6072.8331
90 NA 13309.6531
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 8853.3439 9376.8300
[2,] 11113.0273 8853.3439
[3,] 13026.3051 11113.0273
[4,] 14622.8873 13026.3051
[5,] 13360.1517 14622.8873
[6,] 15837.4748 13360.1517
[7,] 14664.4476 15837.4748
[8,] 12833.6694 14664.4476
[9,] 11045.7397 12833.6694
[10,] 11380.1736 11045.7397
[11,] 8994.1187 11380.1736
[12,] 7094.2574 8994.1187
[13,] 2235.6623 7094.2574
[14,] 388.8190 2235.6623
[15,] -2548.1527 388.8190
[16,] -4217.6259 -2548.1527
[17,] -8284.4447 -4217.6259
[18,] -7232.8722 -8284.4447
[19,] -9085.6499 -7232.8722
[20,] -10722.1230 -9085.6499
[21,] -12498.9695 -10722.1230
[22,] -9948.2029 -12498.9695
[23,] -11972.7865 -9948.2029
[24,] -9357.8162 -11972.7865
[25,] -12261.4686 -9357.8162
[26,] -8575.8406 -12261.4686
[27,] -10008.3965 -8575.8406
[28,] -8693.7311 -10008.3965
[29,] -11632.4667 -8693.7311
[30,] -8133.7555 -11632.4667
[31,] -10328.4223 -8133.7555
[32,] -11633.0618 -10328.4223
[33,] -13674.8529 -11633.0618
[34,] -10895.2249 -13674.8529
[35,] -11244.6976 -10895.2249
[36,] -10623.1985 -11244.6976
[37,] -11265.7678 -10623.1985
[38,] -8453.1121 -11265.7678
[39,] -8759.7788 -8453.1121
[40,] -9536.3075 -8759.7788
[41,] -8978.9045 -9536.3075
[42,] -8266.5814 -8978.9045
[43,] -8688.3591 -8266.5814
[44,] -8427.7769 -8688.3591
[45,] -9398.6511 -8427.7769
[46,] -6784.0508 -9398.6511
[47,] -6727.6622 -6784.0508
[48,] -3834.7472 -6727.6622
[49,] -2312.2887 -3834.7472
[50,] -1752.5776 -2312.2887
[51,] -103.3275 -1752.5776
[52,] 2089.1993 -103.3275
[53,] 646.2973 2089.1993
[54,] 545.9530 646.2973
[55,] -280.7692 545.9530
[56,] 2141.0071 -280.7692
[57,] 2018.3824 2141.0071
[58,] 3181.0936 2018.3824
[59,] 5401.7318 3181.0936
[60,] 6020.4249 5401.7318
[61,] 7659.8002 6020.4249
[62,] 6839.3727 7659.8002
[63,] 7298.5119 6839.3727
[64,] 5681.9555 7298.5119
[65,] 8093.3308 5681.9555
[66,] 7665.9311 8093.3308
[67,] 9439.2089 7665.9311
[68,] 7529.6802 9439.2089
[69,] 9634.9446 7529.6802
[70,] 8236.3785 9634.9446
[71,] 9364.4622 8236.3785
[72,] 7088.0445 9364.4622
[73,] 8585.3921 7088.0445
[74,] 6667.9369 8585.3921
[75,] 9774.2146 6667.9369
[76,] 6465.7969 9774.2146
[77,] 6193.8118 6465.7969
[78,] 3026.0793 6193.8118
[79,] 4754.0799 3026.0793
[80,] 2541.3572 4754.0799
[81,] 2137.3998 2541.3572
[82,] 1650.9169 2137.3998
[83,] 1984.7510 1650.9169
[84,] -365.8885 1984.7510
[85,] -3466.0954 -365.8885
[86,] -6490.9387 -3466.0954
[87,] -4953.8293 -6490.9387
[88,] -6072.8331 -4953.8293
[89,] 13309.6531 -6072.8331
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 8853.3439 9376.8300
2 11113.0273 8853.3439
3 13026.3051 11113.0273
4 14622.8873 13026.3051
5 13360.1517 14622.8873
6 15837.4748 13360.1517
7 14664.4476 15837.4748
8 12833.6694 14664.4476
9 11045.7397 12833.6694
10 11380.1736 11045.7397
11 8994.1187 11380.1736
12 7094.2574 8994.1187
13 2235.6623 7094.2574
14 388.8190 2235.6623
15 -2548.1527 388.8190
16 -4217.6259 -2548.1527
17 -8284.4447 -4217.6259
18 -7232.8722 -8284.4447
19 -9085.6499 -7232.8722
20 -10722.1230 -9085.6499
21 -12498.9695 -10722.1230
22 -9948.2029 -12498.9695
23 -11972.7865 -9948.2029
24 -9357.8162 -11972.7865
25 -12261.4686 -9357.8162
26 -8575.8406 -12261.4686
27 -10008.3965 -8575.8406
28 -8693.7311 -10008.3965
29 -11632.4667 -8693.7311
30 -8133.7555 -11632.4667
31 -10328.4223 -8133.7555
32 -11633.0618 -10328.4223
33 -13674.8529 -11633.0618
34 -10895.2249 -13674.8529
35 -11244.6976 -10895.2249
36 -10623.1985 -11244.6976
37 -11265.7678 -10623.1985
38 -8453.1121 -11265.7678
39 -8759.7788 -8453.1121
40 -9536.3075 -8759.7788
41 -8978.9045 -9536.3075
42 -8266.5814 -8978.9045
43 -8688.3591 -8266.5814
44 -8427.7769 -8688.3591
45 -9398.6511 -8427.7769
46 -6784.0508 -9398.6511
47 -6727.6622 -6784.0508
48 -3834.7472 -6727.6622
49 -2312.2887 -3834.7472
50 -1752.5776 -2312.2887
51 -103.3275 -1752.5776
52 2089.1993 -103.3275
53 646.2973 2089.1993
54 545.9530 646.2973
55 -280.7692 545.9530
56 2141.0071 -280.7692
57 2018.3824 2141.0071
58 3181.0936 2018.3824
59 5401.7318 3181.0936
60 6020.4249 5401.7318
61 7659.8002 6020.4249
62 6839.3727 7659.8002
63 7298.5119 6839.3727
64 5681.9555 7298.5119
65 8093.3308 5681.9555
66 7665.9311 8093.3308
67 9439.2089 7665.9311
68 7529.6802 9439.2089
69 9634.9446 7529.6802
70 8236.3785 9634.9446
71 9364.4622 8236.3785
72 7088.0445 9364.4622
73 8585.3921 7088.0445
74 6667.9369 8585.3921
75 9774.2146 6667.9369
76 6465.7969 9774.2146
77 6193.8118 6465.7969
78 3026.0793 6193.8118
79 4754.0799 3026.0793
80 2541.3572 4754.0799
81 2137.3998 2541.3572
82 1650.9169 2137.3998
83 1984.7510 1650.9169
84 -365.8885 1984.7510
85 -3466.0954 -365.8885
86 -6490.9387 -3466.0954
87 -4953.8293 -6490.9387
88 -6072.8331 -4953.8293
89 13309.6531 -6072.8331
> 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/7aidj1258724451.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/83p5n1258724451.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/9tb311258724451.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/108o3n1258724451.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/1183c51258724451.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/12yen41258724451.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/131wa51258724451.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/14ctm61258724451.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/156ovl1258724451.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/16de381258724451.tab")
+ }
>
> system("convert tmp/191ti1258724451.ps tmp/191ti1258724451.png")
> system("convert tmp/21lkh1258724451.ps tmp/21lkh1258724451.png")
> system("convert tmp/328qj1258724451.ps tmp/328qj1258724451.png")
> system("convert tmp/409dh1258724451.ps tmp/409dh1258724451.png")
> system("convert tmp/5n8bg1258724451.ps tmp/5n8bg1258724451.png")
> system("convert tmp/6buxr1258724451.ps tmp/6buxr1258724451.png")
> system("convert tmp/7aidj1258724451.ps tmp/7aidj1258724451.png")
> system("convert tmp/83p5n1258724451.ps tmp/83p5n1258724451.png")
> system("convert tmp/9tb311258724451.ps tmp/9tb311258724451.png")
> system("convert tmp/108o3n1258724451.ps tmp/108o3n1258724451.png")
>
>
> proc.time()
user system elapsed
2.794 1.617 3.215