R version 2.7.0 (2008-04-22)
Copyright (C) 2008 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> x <- array(list(3353,0,3480,0,3098,0,2944,0,3389,0,3497,0,4404,0,3849,0,3734,0,3060,0,3507,0,3287,0,3215,0,3764,0,2734,0,2837,0,2766,0,3851,0,3289,0,3848,0,3348,0,3682,0,4058,0,3655,1,3811,1,3341,1,3032,1,3475,1,3353,1,3186,1,3902,1,4164,1,3499,1,4145,1,3796,1,3711,1,3949,1,3740,1,3243,1,4407,1,4814,1,3908,1,5250,1,3937,1,4004,1,5560,1,3922,1,3759,1,4138,1,4634,1,3996,1,4308,1,4142,1,4429,1,5219,1,4929,1,5754,1,5592,1,4163,1,4962,1,5208,1,4755,1,4491,1,5732,1,5730,1,5024,1,6056,1,4901,1,5353,1,5578,1,4618,1,4724,1,5011,1,5298,1,4143,1,4617,1,4727,1,4207,1,5112,1,4190,1,4098,1,5071,1,4177,1,4598,1,3757,1,5591,1,4218,1,3780,1,4336,1,4870,1,4422,1,4727,1,4459,1),dim=c(2,93),dimnames=list(c('y','d'),1:93))
> y <- array(NA,dim=c(2,93),dimnames=list(c('y','d'),1:93))
> for (i in 1:dim(x)[1])
+ {
+ for (j in 1:dim(x)[2])
+ {
+ y[i,j] <- as.numeric(x[i,j])
+ }
+ }
> par3 = 'Linear Trend'
> par2 = 'Include Monthly 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 d M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 3353 0 1 0 0 0 0 0 0 0 0 0 0 1
2 3480 0 0 1 0 0 0 0 0 0 0 0 0 2
3 3098 0 0 0 1 0 0 0 0 0 0 0 0 3
4 2944 0 0 0 0 1 0 0 0 0 0 0 0 4
5 3389 0 0 0 0 0 1 0 0 0 0 0 0 5
6 3497 0 0 0 0 0 0 1 0 0 0 0 0 6
7 4404 0 0 0 0 0 0 0 1 0 0 0 0 7
8 3849 0 0 0 0 0 0 0 0 1 0 0 0 8
9 3734 0 0 0 0 0 0 0 0 0 1 0 0 9
10 3060 0 0 0 0 0 0 0 0 0 0 1 0 10
11 3507 0 0 0 0 0 0 0 0 0 0 0 1 11
12 3287 0 0 0 0 0 0 0 0 0 0 0 0 12
13 3215 0 1 0 0 0 0 0 0 0 0 0 0 13
14 3764 0 0 1 0 0 0 0 0 0 0 0 0 14
15 2734 0 0 0 1 0 0 0 0 0 0 0 0 15
16 2837 0 0 0 0 1 0 0 0 0 0 0 0 16
17 2766 0 0 0 0 0 1 0 0 0 0 0 0 17
18 3851 0 0 0 0 0 0 1 0 0 0 0 0 18
19 3289 0 0 0 0 0 0 0 1 0 0 0 0 19
20 3848 0 0 0 0 0 0 0 0 1 0 0 0 20
21 3348 0 0 0 0 0 0 0 0 0 1 0 0 21
22 3682 0 0 0 0 0 0 0 0 0 0 1 0 22
23 4058 0 0 0 0 0 0 0 0 0 0 0 1 23
24 3655 1 0 0 0 0 0 0 0 0 0 0 0 24
25 3811 1 1 0 0 0 0 0 0 0 0 0 0 25
26 3341 1 0 1 0 0 0 0 0 0 0 0 0 26
27 3032 1 0 0 1 0 0 0 0 0 0 0 0 27
28 3475 1 0 0 0 1 0 0 0 0 0 0 0 28
29 3353 1 0 0 0 0 1 0 0 0 0 0 0 29
30 3186 1 0 0 0 0 0 1 0 0 0 0 0 30
31 3902 1 0 0 0 0 0 0 1 0 0 0 0 31
32 4164 1 0 0 0 0 0 0 0 1 0 0 0 32
33 3499 1 0 0 0 0 0 0 0 0 1 0 0 33
34 4145 1 0 0 0 0 0 0 0 0 0 1 0 34
35 3796 1 0 0 0 0 0 0 0 0 0 0 1 35
36 3711 1 0 0 0 0 0 0 0 0 0 0 0 36
37 3949 1 1 0 0 0 0 0 0 0 0 0 0 37
38 3740 1 0 1 0 0 0 0 0 0 0 0 0 38
39 3243 1 0 0 1 0 0 0 0 0 0 0 0 39
40 4407 1 0 0 0 1 0 0 0 0 0 0 0 40
41 4814 1 0 0 0 0 1 0 0 0 0 0 0 41
42 3908 1 0 0 0 0 0 1 0 0 0 0 0 42
43 5250 1 0 0 0 0 0 0 1 0 0 0 0 43
44 3937 1 0 0 0 0 0 0 0 1 0 0 0 44
45 4004 1 0 0 0 0 0 0 0 0 1 0 0 45
46 5560 1 0 0 0 0 0 0 0 0 0 1 0 46
47 3922 1 0 0 0 0 0 0 0 0 0 0 1 47
48 3759 1 0 0 0 0 0 0 0 0 0 0 0 48
49 4138 1 1 0 0 0 0 0 0 0 0 0 0 49
50 4634 1 0 1 0 0 0 0 0 0 0 0 0 50
51 3996 1 0 0 1 0 0 0 0 0 0 0 0 51
52 4308 1 0 0 0 1 0 0 0 0 0 0 0 52
53 4142 1 0 0 0 0 1 0 0 0 0 0 0 53
54 4429 1 0 0 0 0 0 1 0 0 0 0 0 54
55 5219 1 0 0 0 0 0 0 1 0 0 0 0 55
56 4929 1 0 0 0 0 0 0 0 1 0 0 0 56
57 5754 1 0 0 0 0 0 0 0 0 1 0 0 57
58 5592 1 0 0 0 0 0 0 0 0 0 1 0 58
59 4163 1 0 0 0 0 0 0 0 0 0 0 1 59
60 4962 1 0 0 0 0 0 0 0 0 0 0 0 60
61 5208 1 1 0 0 0 0 0 0 0 0 0 0 61
62 4755 1 0 1 0 0 0 0 0 0 0 0 0 62
63 4491 1 0 0 1 0 0 0 0 0 0 0 0 63
64 5732 1 0 0 0 1 0 0 0 0 0 0 0 64
65 5730 1 0 0 0 0 1 0 0 0 0 0 0 65
66 5024 1 0 0 0 0 0 1 0 0 0 0 0 66
67 6056 1 0 0 0 0 0 0 1 0 0 0 0 67
68 4901 1 0 0 0 0 0 0 0 1 0 0 0 68
69 5353 1 0 0 0 0 0 0 0 0 1 0 0 69
70 5578 1 0 0 0 0 0 0 0 0 0 1 0 70
71 4618 1 0 0 0 0 0 0 0 0 0 0 1 71
72 4724 1 0 0 0 0 0 0 0 0 0 0 0 72
73 5011 1 1 0 0 0 0 0 0 0 0 0 0 73
74 5298 1 0 1 0 0 0 0 0 0 0 0 0 74
75 4143 1 0 0 1 0 0 0 0 0 0 0 0 75
76 4617 1 0 0 0 1 0 0 0 0 0 0 0 76
77 4727 1 0 0 0 0 1 0 0 0 0 0 0 77
78 4207 1 0 0 0 0 0 1 0 0 0 0 0 78
79 5112 1 0 0 0 0 0 0 1 0 0 0 0 79
80 4190 1 0 0 0 0 0 0 0 1 0 0 0 80
81 4098 1 0 0 0 0 0 0 0 0 1 0 0 81
82 5071 1 0 0 0 0 0 0 0 0 0 1 0 82
83 4177 1 0 0 0 0 0 0 0 0 0 0 1 83
84 4598 1 0 0 0 0 0 0 0 0 0 0 0 84
85 3757 1 1 0 0 0 0 0 0 0 0 0 0 85
86 5591 1 0 1 0 0 0 0 0 0 0 0 0 86
87 4218 1 0 0 1 0 0 0 0 0 0 0 0 87
88 3780 1 0 0 0 1 0 0 0 0 0 0 0 88
89 4336 1 0 0 0 0 1 0 0 0 0 0 0 89
90 4870 1 0 0 0 0 0 1 0 0 0 0 0 90
91 4422 1 0 0 0 0 0 0 1 0 0 0 0 91
92 4727 1 0 0 0 0 0 0 0 1 0 0 0 92
93 4459 1 0 0 0 0 0 0 0 0 1 0 0 93
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) d M1 M2 M3 M4
3094.418 305.913 65.973 320.623 -400.852 -23.202
M5 M6 M7 M8 M9 M10
105.948 54.848 624.623 220.523 168.048 644.938
M11 t
-5.823 15.475
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-1024.68 -364.75 -17.10 358.01 1364.47
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3094.418 240.054 12.891 < 2e-16 ***
d 305.913 202.338 1.512 0.1346
M1 65.973 287.102 0.230 0.8189
M2 320.623 287.121 1.117 0.2675
M3 -400.852 287.176 -1.396 0.1667
M4 -23.202 287.269 -0.081 0.9358
M5 105.948 287.398 0.369 0.7134
M6 54.848 287.564 0.191 0.8492
M7 624.623 287.767 2.171 0.0330 *
M8 220.523 288.006 0.766 0.4461
M9 168.048 288.282 0.583 0.5616
M10 644.938 297.153 2.170 0.0330 *
M11 -5.823 297.337 -0.020 0.9844
t 15.475 3.255 4.755 8.79e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 554 on 79 degrees of freedom
Multiple R-squared: 0.5728, Adjusted R-squared: 0.5025
F-statistic: 8.148 on 13 and 79 DF, p-value: 4.3e-10
> 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.116397292 0.232794585 0.8836027
[2,] 0.092545090 0.185090180 0.9074549
[3,] 0.194856577 0.389713155 0.8051434
[4,] 0.117412243 0.234824487 0.8825878
[5,] 0.065136609 0.130273218 0.9348634
[6,] 0.095155394 0.190310788 0.9048446
[7,] 0.090207687 0.180415375 0.9097923
[8,] 0.052278434 0.104556867 0.9477216
[9,] 0.029738195 0.059476389 0.9702618
[10,] 0.033090831 0.066181661 0.9669092
[11,] 0.019079810 0.038159619 0.9809202
[12,] 0.014295071 0.028590141 0.9857049
[13,] 0.009249882 0.018499763 0.9907501
[14,] 0.014232712 0.028465424 0.9857673
[15,] 0.010037977 0.020075955 0.9899620
[16,] 0.005759788 0.011519576 0.9942402
[17,] 0.004255988 0.008511975 0.9957440
[18,] 0.006781782 0.013563565 0.9932182
[19,] 0.004007608 0.008015217 0.9959924
[20,] 0.002710348 0.005420697 0.9972897
[21,] 0.002220037 0.004440073 0.9977800
[22,] 0.002278909 0.004557817 0.9977211
[23,] 0.001913278 0.003826557 0.9980867
[24,] 0.010570361 0.021140722 0.9894296
[25,] 0.047915518 0.095831037 0.9520845
[26,] 0.039358515 0.078717029 0.9606415
[27,] 0.057176156 0.114352311 0.9428238
[28,] 0.059248567 0.118497135 0.9407514
[29,] 0.063314399 0.126628798 0.9366856
[30,] 0.154358581 0.308717162 0.8456414
[31,] 0.136846304 0.273692608 0.8631537
[32,] 0.173562032 0.347124063 0.8264380
[33,] 0.164065113 0.328130227 0.8359349
[34,] 0.182294038 0.364588075 0.8177060
[35,] 0.171092580 0.342185160 0.8289074
[36,] 0.174285275 0.348570551 0.8257147
[37,] 0.271475448 0.542950896 0.7285246
[38,] 0.297119023 0.594238047 0.7028810
[39,] 0.292702099 0.585404198 0.7072979
[40,] 0.254948588 0.509897176 0.7450514
[41,] 0.373884257 0.747768513 0.6261157
[42,] 0.345126200 0.690252401 0.6548738
[43,] 0.376300968 0.752601937 0.6236990
[44,] 0.330566002 0.661132005 0.6694340
[45,] 0.281980800 0.563961599 0.7180192
[46,] 0.479573823 0.959147645 0.5204262
[47,] 0.428410733 0.856821466 0.5715893
[48,] 0.563670794 0.872658413 0.4363292
[49,] 0.599723550 0.800552900 0.4002764
[50,] 0.513055227 0.973889547 0.4869448
[51,] 0.572776005 0.854447989 0.4272240
[52,] 0.495049273 0.990098546 0.5049507
[53,] 0.524482831 0.951034339 0.4755172
[54,] 0.444725829 0.889451658 0.5552742
[55,] 0.372151620 0.744303240 0.6278484
[56,] 0.281336111 0.562672222 0.7186639
[57,] 0.484577014 0.969154028 0.5154230
[58,] 0.381250610 0.762501219 0.6187494
[59,] 0.273681611 0.547363222 0.7263184
[60,] 0.335950942 0.671901884 0.6640491
> postscript(file="/var/www/html/rcomp/tmp/1xftf1229180879.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/2p73s1229180879.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/397081229180879.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/4u69j1229180879.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/5y67e1229180879.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 = 93
Frequency = 1
1 2 3 4 5 6
177.134238 34.009238 358.009238 -189.115762 111.259238 254.884238
7 8 9 10 11 12
576.634238 410.259238 332.259238 -834.105431 248.180283 6.882155
13 14 15 16 17 18
-146.565593 132.309407 -191.690593 -481.815593 -697.440593 423.184407
19 20 21 22 23 24
-724.065593 223.559407 -239.440593 -397.805262 613.480452 -116.730782
25 26 27 28 29 30
-42.178530 -782.303530 -385.303530 -335.428530 -602.053530 -733.428530
31 32 33 34 35 36
-602.678530 47.946470 -580.053530 -426.418199 -140.132485 -246.430613
37 38 39 40 41 42
-89.878361 -569.003361 -360.003361 410.871639 673.246639 -197.128361
43 44 45 46 47 48
559.621639 -364.753361 -260.753361 802.881970 -199.832316 -384.130444
49 50 51 52 53 54
-86.578192 139.296808 207.296808 126.171808 -184.453192 138.171808
55 56 57 58 59 60
342.921808 441.546808 1303.546808 649.182139 -144.532147 633.169725
61 62 63 64 65 66
797.721977 74.596977 516.596977 1364.471977 1217.846977 547.471977
67 68 69 70 71 72
994.221977 227.846977 716.846977 449.482308 124.768022 209.469894
73 74 75 76 77 78
415.022146 431.897146 -17.102854 63.772146 29.147146 -455.227854
79 80 81 82 83 84
-135.477854 -668.852854 -723.852854 -243.217523 -501.931809 -102.229936
85 86 87 88 89 90
-1024.677685 539.197315 -127.802685 -958.927685 -547.552685 22.072315
91 92 93
-1011.177685 -317.552685 -548.552685
> postscript(file="/var/www/html/rcomp/tmp/695091229180879.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 = 93
Frequency = 1
lag(myerror, k = 1) myerror
0 177.134238 NA
1 34.009238 177.134238
2 358.009238 34.009238
3 -189.115762 358.009238
4 111.259238 -189.115762
5 254.884238 111.259238
6 576.634238 254.884238
7 410.259238 576.634238
8 332.259238 410.259238
9 -834.105431 332.259238
10 248.180283 -834.105431
11 6.882155 248.180283
12 -146.565593 6.882155
13 132.309407 -146.565593
14 -191.690593 132.309407
15 -481.815593 -191.690593
16 -697.440593 -481.815593
17 423.184407 -697.440593
18 -724.065593 423.184407
19 223.559407 -724.065593
20 -239.440593 223.559407
21 -397.805262 -239.440593
22 613.480452 -397.805262
23 -116.730782 613.480452
24 -42.178530 -116.730782
25 -782.303530 -42.178530
26 -385.303530 -782.303530
27 -335.428530 -385.303530
28 -602.053530 -335.428530
29 -733.428530 -602.053530
30 -602.678530 -733.428530
31 47.946470 -602.678530
32 -580.053530 47.946470
33 -426.418199 -580.053530
34 -140.132485 -426.418199
35 -246.430613 -140.132485
36 -89.878361 -246.430613
37 -569.003361 -89.878361
38 -360.003361 -569.003361
39 410.871639 -360.003361
40 673.246639 410.871639
41 -197.128361 673.246639
42 559.621639 -197.128361
43 -364.753361 559.621639
44 -260.753361 -364.753361
45 802.881970 -260.753361
46 -199.832316 802.881970
47 -384.130444 -199.832316
48 -86.578192 -384.130444
49 139.296808 -86.578192
50 207.296808 139.296808
51 126.171808 207.296808
52 -184.453192 126.171808
53 138.171808 -184.453192
54 342.921808 138.171808
55 441.546808 342.921808
56 1303.546808 441.546808
57 649.182139 1303.546808
58 -144.532147 649.182139
59 633.169725 -144.532147
60 797.721977 633.169725
61 74.596977 797.721977
62 516.596977 74.596977
63 1364.471977 516.596977
64 1217.846977 1364.471977
65 547.471977 1217.846977
66 994.221977 547.471977
67 227.846977 994.221977
68 716.846977 227.846977
69 449.482308 716.846977
70 124.768022 449.482308
71 209.469894 124.768022
72 415.022146 209.469894
73 431.897146 415.022146
74 -17.102854 431.897146
75 63.772146 -17.102854
76 29.147146 63.772146
77 -455.227854 29.147146
78 -135.477854 -455.227854
79 -668.852854 -135.477854
80 -723.852854 -668.852854
81 -243.217523 -723.852854
82 -501.931809 -243.217523
83 -102.229936 -501.931809
84 -1024.677685 -102.229936
85 539.197315 -1024.677685
86 -127.802685 539.197315
87 -958.927685 -127.802685
88 -547.552685 -958.927685
89 22.072315 -547.552685
90 -1011.177685 22.072315
91 -317.552685 -1011.177685
92 -548.552685 -317.552685
93 NA -548.552685
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 34.009238 177.134238
[2,] 358.009238 34.009238
[3,] -189.115762 358.009238
[4,] 111.259238 -189.115762
[5,] 254.884238 111.259238
[6,] 576.634238 254.884238
[7,] 410.259238 576.634238
[8,] 332.259238 410.259238
[9,] -834.105431 332.259238
[10,] 248.180283 -834.105431
[11,] 6.882155 248.180283
[12,] -146.565593 6.882155
[13,] 132.309407 -146.565593
[14,] -191.690593 132.309407
[15,] -481.815593 -191.690593
[16,] -697.440593 -481.815593
[17,] 423.184407 -697.440593
[18,] -724.065593 423.184407
[19,] 223.559407 -724.065593
[20,] -239.440593 223.559407
[21,] -397.805262 -239.440593
[22,] 613.480452 -397.805262
[23,] -116.730782 613.480452
[24,] -42.178530 -116.730782
[25,] -782.303530 -42.178530
[26,] -385.303530 -782.303530
[27,] -335.428530 -385.303530
[28,] -602.053530 -335.428530
[29,] -733.428530 -602.053530
[30,] -602.678530 -733.428530
[31,] 47.946470 -602.678530
[32,] -580.053530 47.946470
[33,] -426.418199 -580.053530
[34,] -140.132485 -426.418199
[35,] -246.430613 -140.132485
[36,] -89.878361 -246.430613
[37,] -569.003361 -89.878361
[38,] -360.003361 -569.003361
[39,] 410.871639 -360.003361
[40,] 673.246639 410.871639
[41,] -197.128361 673.246639
[42,] 559.621639 -197.128361
[43,] -364.753361 559.621639
[44,] -260.753361 -364.753361
[45,] 802.881970 -260.753361
[46,] -199.832316 802.881970
[47,] -384.130444 -199.832316
[48,] -86.578192 -384.130444
[49,] 139.296808 -86.578192
[50,] 207.296808 139.296808
[51,] 126.171808 207.296808
[52,] -184.453192 126.171808
[53,] 138.171808 -184.453192
[54,] 342.921808 138.171808
[55,] 441.546808 342.921808
[56,] 1303.546808 441.546808
[57,] 649.182139 1303.546808
[58,] -144.532147 649.182139
[59,] 633.169725 -144.532147
[60,] 797.721977 633.169725
[61,] 74.596977 797.721977
[62,] 516.596977 74.596977
[63,] 1364.471977 516.596977
[64,] 1217.846977 1364.471977
[65,] 547.471977 1217.846977
[66,] 994.221977 547.471977
[67,] 227.846977 994.221977
[68,] 716.846977 227.846977
[69,] 449.482308 716.846977
[70,] 124.768022 449.482308
[71,] 209.469894 124.768022
[72,] 415.022146 209.469894
[73,] 431.897146 415.022146
[74,] -17.102854 431.897146
[75,] 63.772146 -17.102854
[76,] 29.147146 63.772146
[77,] -455.227854 29.147146
[78,] -135.477854 -455.227854
[79,] -668.852854 -135.477854
[80,] -723.852854 -668.852854
[81,] -243.217523 -723.852854
[82,] -501.931809 -243.217523
[83,] -102.229936 -501.931809
[84,] -1024.677685 -102.229936
[85,] 539.197315 -1024.677685
[86,] -127.802685 539.197315
[87,] -958.927685 -127.802685
[88,] -547.552685 -958.927685
[89,] 22.072315 -547.552685
[90,] -1011.177685 22.072315
[91,] -317.552685 -1011.177685
[92,] -548.552685 -317.552685
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 34.009238 177.134238
2 358.009238 34.009238
3 -189.115762 358.009238
4 111.259238 -189.115762
5 254.884238 111.259238
6 576.634238 254.884238
7 410.259238 576.634238
8 332.259238 410.259238
9 -834.105431 332.259238
10 248.180283 -834.105431
11 6.882155 248.180283
12 -146.565593 6.882155
13 132.309407 -146.565593
14 -191.690593 132.309407
15 -481.815593 -191.690593
16 -697.440593 -481.815593
17 423.184407 -697.440593
18 -724.065593 423.184407
19 223.559407 -724.065593
20 -239.440593 223.559407
21 -397.805262 -239.440593
22 613.480452 -397.805262
23 -116.730782 613.480452
24 -42.178530 -116.730782
25 -782.303530 -42.178530
26 -385.303530 -782.303530
27 -335.428530 -385.303530
28 -602.053530 -335.428530
29 -733.428530 -602.053530
30 -602.678530 -733.428530
31 47.946470 -602.678530
32 -580.053530 47.946470
33 -426.418199 -580.053530
34 -140.132485 -426.418199
35 -246.430613 -140.132485
36 -89.878361 -246.430613
37 -569.003361 -89.878361
38 -360.003361 -569.003361
39 410.871639 -360.003361
40 673.246639 410.871639
41 -197.128361 673.246639
42 559.621639 -197.128361
43 -364.753361 559.621639
44 -260.753361 -364.753361
45 802.881970 -260.753361
46 -199.832316 802.881970
47 -384.130444 -199.832316
48 -86.578192 -384.130444
49 139.296808 -86.578192
50 207.296808 139.296808
51 126.171808 207.296808
52 -184.453192 126.171808
53 138.171808 -184.453192
54 342.921808 138.171808
55 441.546808 342.921808
56 1303.546808 441.546808
57 649.182139 1303.546808
58 -144.532147 649.182139
59 633.169725 -144.532147
60 797.721977 633.169725
61 74.596977 797.721977
62 516.596977 74.596977
63 1364.471977 516.596977
64 1217.846977 1364.471977
65 547.471977 1217.846977
66 994.221977 547.471977
67 227.846977 994.221977
68 716.846977 227.846977
69 449.482308 716.846977
70 124.768022 449.482308
71 209.469894 124.768022
72 415.022146 209.469894
73 431.897146 415.022146
74 -17.102854 431.897146
75 63.772146 -17.102854
76 29.147146 63.772146
77 -455.227854 29.147146
78 -135.477854 -455.227854
79 -668.852854 -135.477854
80 -723.852854 -668.852854
81 -243.217523 -723.852854
82 -501.931809 -243.217523
83 -102.229936 -501.931809
84 -1024.677685 -102.229936
85 539.197315 -1024.677685
86 -127.802685 539.197315
87 -958.927685 -127.802685
88 -547.552685 -958.927685
89 22.072315 -547.552685
90 -1011.177685 22.072315
91 -317.552685 -1011.177685
92 -548.552685 -317.552685
> 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/74uo51229180879.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/8ib7b1229180879.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/9wzwy1229180879.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/10iev41229180879.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/11t1fn1229180879.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/12jruz1229180880.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/13hwpt1229180880.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/14ks8g1229180880.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/15or501229180880.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/16peeo1229180880.tab")
+ }
>
> system("convert tmp/1xftf1229180879.ps tmp/1xftf1229180879.png")
> system("convert tmp/2p73s1229180879.ps tmp/2p73s1229180879.png")
> system("convert tmp/397081229180879.ps tmp/397081229180879.png")
> system("convert tmp/4u69j1229180879.ps tmp/4u69j1229180879.png")
> system("convert tmp/5y67e1229180879.ps tmp/5y67e1229180879.png")
> system("convert tmp/695091229180879.ps tmp/695091229180879.png")
> system("convert tmp/74uo51229180879.ps tmp/74uo51229180879.png")
> system("convert tmp/8ib7b1229180879.ps tmp/8ib7b1229180879.png")
> system("convert tmp/9wzwy1229180879.ps tmp/9wzwy1229180879.png")
> system("convert tmp/10iev41229180879.ps tmp/10iev41229180879.png")
>
>
> proc.time()
user system elapsed
5.732 2.825 6.123