R version 2.15.2 (2012-10-26) -- "Trick or Treat"
Copyright (C) 2012 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: i686-pc-linux-gnu (32-bit)
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(-3
+ ,14
+ ,24
+ ,6
+ ,17
+ ,-4
+ ,16
+ ,24
+ ,6
+ ,13
+ ,-7
+ ,19
+ ,31
+ ,5
+ ,12
+ ,-7
+ ,18
+ ,25
+ ,5
+ ,13
+ ,-7
+ ,19
+ ,28
+ ,3
+ ,10
+ ,-3
+ ,20
+ ,24
+ ,5
+ ,14
+ ,0
+ ,20
+ ,25
+ ,5
+ ,13
+ ,-5
+ ,24
+ ,16
+ ,5
+ ,10
+ ,-3
+ ,18
+ ,17
+ ,3
+ ,11
+ ,3
+ ,15
+ ,11
+ ,6
+ ,12
+ ,2
+ ,25
+ ,12
+ ,6
+ ,7
+ ,-7
+ ,23
+ ,39
+ ,4
+ ,11
+ ,-1
+ ,20
+ ,19
+ ,6
+ ,9
+ ,0
+ ,20
+ ,14
+ ,5
+ ,13
+ ,-3
+ ,22
+ ,15
+ ,4
+ ,12
+ ,4
+ ,25
+ ,7
+ ,5
+ ,5
+ ,2
+ ,22
+ ,12
+ ,5
+ ,13
+ ,3
+ ,26
+ ,12
+ ,4
+ ,11
+ ,0
+ ,27
+ ,14
+ ,3
+ ,8
+ ,-10
+ ,41
+ ,9
+ ,2
+ ,8
+ ,-10
+ ,29
+ ,8
+ ,3
+ ,8
+ ,-9
+ ,33
+ ,4
+ ,2
+ ,8
+ ,-22
+ ,39
+ ,7
+ ,-1
+ ,0
+ ,-16
+ ,27
+ ,3
+ ,0
+ ,3
+ ,-18
+ ,27
+ ,5
+ ,-2
+ ,0
+ ,-14
+ ,25
+ ,0
+ ,1
+ ,-1
+ ,-12
+ ,19
+ ,-2
+ ,-2
+ ,-1
+ ,-17
+ ,15
+ ,6
+ ,-2
+ ,-4
+ ,-23
+ ,19
+ ,11
+ ,-2
+ ,1
+ ,-28
+ ,23
+ ,9
+ ,-6
+ ,-1
+ ,-31
+ ,23
+ ,17
+ ,-4
+ ,0
+ ,-21
+ ,7
+ ,21
+ ,-2
+ ,-1
+ ,-19
+ ,1
+ ,21
+ ,0
+ ,6
+ ,-22
+ ,7
+ ,41
+ ,-5
+ ,0
+ ,-22
+ ,4
+ ,57
+ ,-4
+ ,-3
+ ,-25
+ ,-8
+ ,65
+ ,-5
+ ,-3
+ ,-16
+ ,-14
+ ,68
+ ,-1
+ ,4
+ ,-22
+ ,-10
+ ,73
+ ,-2
+ ,1
+ ,-21
+ ,-11
+ ,71
+ ,-4
+ ,0
+ ,-10
+ ,-10
+ ,71
+ ,-1
+ ,-4
+ ,-7
+ ,-8
+ ,70
+ ,1
+ ,-2
+ ,-5
+ ,-8
+ ,69
+ ,1
+ ,3
+ ,-4
+ ,-7
+ ,65
+ ,-2
+ ,2
+ ,7
+ ,-8
+ ,57
+ ,1
+ ,5
+ ,6
+ ,-4
+ ,57
+ ,1
+ ,6
+ ,3
+ ,3
+ ,57
+ ,3
+ ,6
+ ,10
+ ,-5
+ ,55
+ ,3
+ ,3
+ ,0
+ ,-4
+ ,65
+ ,1
+ ,4
+ ,-2
+ ,5
+ ,65
+ ,1
+ ,7
+ ,-1
+ ,3
+ ,64
+ ,0
+ ,5
+ ,2
+ ,6
+ ,60
+ ,2
+ ,6
+ ,8
+ ,10
+ ,43
+ ,2
+ ,1
+ ,-6
+ ,16
+ ,47
+ ,-1
+ ,3
+ ,-4
+ ,11
+ ,40
+ ,1
+ ,6
+ ,4
+ ,10
+ ,31
+ ,0
+ ,0
+ ,7
+ ,21
+ ,27
+ ,1
+ ,3
+ ,3
+ ,18
+ ,24
+ ,1
+ ,4
+ ,3
+ ,20
+ ,23
+ ,3
+ ,7
+ ,8
+ ,18
+ ,17
+ ,2
+ ,6
+ ,3
+ ,23
+ ,16
+ ,0
+ ,6
+ ,-3
+ ,28
+ ,15
+ ,0
+ ,6
+ ,4
+ ,31
+ ,8
+ ,3
+ ,6
+ ,-5
+ ,38
+ ,5
+ ,-2
+ ,2
+ ,-1
+ ,27
+ ,6
+ ,0
+ ,2
+ ,5
+ ,21
+ ,5
+ ,1
+ ,2
+ ,0
+ ,31
+ ,12
+ ,-1
+ ,3
+ ,-6
+ ,31
+ ,8
+ ,-2
+ ,-1
+ ,-13
+ ,29
+ ,17
+ ,-1
+ ,-4
+ ,-15
+ ,24
+ ,22
+ ,-1
+ ,4
+ ,-8
+ ,27
+ ,24
+ ,1
+ ,5
+ ,-20
+ ,36
+ ,36
+ ,-2
+ ,3
+ ,-10
+ ,35
+ ,31
+ ,-5
+ ,-1
+ ,-22
+ ,44
+ ,34
+ ,-5
+ ,-4
+ ,-25
+ ,39
+ ,47
+ ,-6
+ ,0
+ ,-10
+ ,26
+ ,33
+ ,-4
+ ,-1
+ ,-8
+ ,27
+ ,35
+ ,-3
+ ,-1
+ ,-9
+ ,17
+ ,31
+ ,-3
+ ,3
+ ,-5
+ ,20
+ ,35
+ ,-1
+ ,2
+ ,-7
+ ,22
+ ,39
+ ,-2
+ ,-4
+ ,-11
+ ,32
+ ,46
+ ,-3
+ ,-3
+ ,-11
+ ,28
+ ,40
+ ,-3
+ ,-1
+ ,-16
+ ,30
+ ,50
+ ,-3
+ ,3)
+ ,dim=c(5
+ ,82)
+ ,dimnames=list(c('Y_t'
+ ,'X_1t'
+ ,'X_2t'
+ ,'X_3t'
+ ,'X_4t')
+ ,1:82))
> y <- array(NA,dim=c(5,82),dimnames=list(c('Y_t','X_1t','X_2t','X_3t','X_4t'),1:82))
> 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'
> 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, 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_t X_1t X_2t X_3t X_4t
1 -3 14 24 6 17
2 -4 16 24 6 13
3 -7 19 31 5 12
4 -7 18 25 5 13
5 -7 19 28 3 10
6 -3 20 24 5 14
7 0 20 25 5 13
8 -5 24 16 5 10
9 -3 18 17 3 11
10 3 15 11 6 12
11 2 25 12 6 7
12 -7 23 39 4 11
13 -1 20 19 6 9
14 0 20 14 5 13
15 -3 22 15 4 12
16 4 25 7 5 5
17 2 22 12 5 13
18 3 26 12 4 11
19 0 27 14 3 8
20 -10 41 9 2 8
21 -10 29 8 3 8
22 -9 33 4 2 8
23 -22 39 7 -1 0
24 -16 27 3 0 3
25 -18 27 5 -2 0
26 -14 25 0 1 -1
27 -12 19 -2 -2 -1
28 -17 15 6 -2 -4
29 -23 19 11 -2 1
30 -28 23 9 -6 -1
31 -31 23 17 -4 0
32 -21 7 21 -2 -1
33 -19 1 21 0 6
34 -22 7 41 -5 0
35 -22 4 57 -4 -3
36 -25 -8 65 -5 -3
37 -16 -14 68 -1 4
38 -22 -10 73 -2 1
39 -21 -11 71 -4 0
40 -10 -10 71 -1 -4
41 -7 -8 70 1 -2
42 -5 -8 69 1 3
43 -4 -7 65 -2 2
44 7 -8 57 1 5
45 6 -4 57 1 6
46 3 3 57 3 6
47 10 -5 55 3 3
48 0 -4 65 1 4
49 -2 5 65 1 7
50 -1 3 64 0 5
51 2 6 60 2 6
52 8 10 43 2 1
53 -6 16 47 -1 3
54 -4 11 40 1 6
55 4 10 31 0 0
56 7 21 27 1 3
57 3 18 24 1 4
58 3 20 23 3 7
59 8 18 17 2 6
60 3 23 16 0 6
61 -3 28 15 0 6
62 4 31 8 3 6
63 -5 38 5 -2 2
64 -1 27 6 0 2
65 5 21 5 1 2
66 0 31 12 -1 3
67 -6 31 8 -2 -1
68 -13 29 17 -1 -4
69 -15 24 22 -1 4
70 -8 27 24 1 5
71 -20 36 36 -2 3
72 -10 35 31 -5 -1
73 -22 44 34 -5 -4
74 -25 39 47 -6 0
75 -10 26 33 -4 -1
76 -8 27 35 -3 -1
77 -9 17 31 -3 3
78 -5 20 35 -1 2
79 -7 22 39 -2 -4
80 -11 32 46 -3 -3
81 -11 28 40 -3 -1
82 -16 30 50 -3 3
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X_1t X_2t X_3t X_4t
-9.9854 0.0673 0.0945 2.9820 -0.5952
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-12.3189 -5.8936 -0.5269 5.5862 13.4964
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -9.98541 3.23389 -3.088 0.0028 **
X_1t 0.06730 0.08732 0.771 0.4432
X_2t 0.09450 0.05830 1.621 0.1091
X_3t 2.98205 0.47531 6.274 1.9e-08 ***
X_4t -0.59517 0.28747 -2.070 0.0418 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 7.018 on 77 degrees of freedom
Multiple R-squared: 0.5028, Adjusted R-squared: 0.477
F-statistic: 19.47 on 4 and 77 DF, p-value: 4.218e-11
> 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,] 8.332781e-02 1.666556e-01 0.9166721918
[2,] 3.138292e-02 6.276585e-02 0.9686170761
[3,] 2.248295e-02 4.496590e-02 0.9775170497
[4,] 1.138240e-02 2.276479e-02 0.9886176034
[5,] 6.754276e-03 1.350855e-02 0.9932457242
[6,] 2.678662e-03 5.357324e-03 0.9973213379
[7,] 9.612553e-04 1.922511e-03 0.9990387447
[8,] 3.518399e-04 7.036798e-04 0.9996481601
[9,] 1.584843e-04 3.169686e-04 0.9998415157
[10,] 7.026007e-05 1.405201e-04 0.9999297399
[11,] 3.907486e-05 7.814971e-05 0.9999609251
[12,] 1.231114e-05 2.462228e-05 0.9999876889
[13,] 1.775694e-04 3.551387e-04 0.9998224306
[14,] 1.505705e-03 3.011410e-03 0.9984942952
[15,] 1.604818e-03 3.209636e-03 0.9983951818
[16,] 6.403441e-03 1.280688e-02 0.9935965593
[17,] 6.049206e-03 1.209841e-02 0.9939507941
[18,] 3.808640e-03 7.617281e-03 0.9961913596
[19,] 4.738369e-03 9.476738e-03 0.9952616308
[20,] 4.533110e-03 9.066220e-03 0.9954668899
[21,] 2.662874e-03 5.325747e-03 0.9973371265
[22,] 3.559227e-03 7.118453e-03 0.9964407734
[23,] 2.088506e-03 4.177013e-03 0.9979114937
[24,] 3.975763e-03 7.951526e-03 0.9960242372
[25,] 4.025642e-03 8.051284e-03 0.9959743580
[26,] 1.630817e-02 3.261633e-02 0.9836918343
[27,] 5.751060e-02 1.150212e-01 0.9424894029
[28,] 6.357275e-02 1.271455e-01 0.9364272459
[29,] 5.662807e-02 1.132561e-01 0.9433719285
[30,] 6.637641e-02 1.327528e-01 0.9336235870
[31,] 1.317643e-01 2.635286e-01 0.8682357010
[32,] 2.197786e-01 4.395572e-01 0.7802213795
[33,] 2.756760e-01 5.513520e-01 0.7243239964
[34,] 3.520655e-01 7.041310e-01 0.6479344956
[35,] 4.631614e-01 9.263228e-01 0.5368385945
[36,] 7.217065e-01 5.565870e-01 0.2782934985
[37,] 8.852610e-01 2.294780e-01 0.1147389794
[38,] 9.406987e-01 1.186026e-01 0.0593013154
[39,] 9.296626e-01 1.406747e-01 0.0703373651
[40,] 9.355853e-01 1.288293e-01 0.0644146699
[41,] 9.277356e-01 1.445289e-01 0.0722644379
[42,] 9.115948e-01 1.768104e-01 0.0884051907
[43,] 9.043615e-01 1.912769e-01 0.0956384692
[44,] 8.827638e-01 2.344723e-01 0.1172361554
[45,] 9.111733e-01 1.776534e-01 0.0888267125
[46,] 8.934242e-01 2.131516e-01 0.1065758014
[47,] 8.869272e-01 2.261457e-01 0.1130728392
[48,] 9.026591e-01 1.946819e-01 0.0973409411
[49,] 9.580249e-01 8.395025e-02 0.0419751268
[50,] 9.562984e-01 8.740317e-02 0.0437015848
[51,] 9.399878e-01 1.200244e-01 0.0600122015
[52,] 9.572495e-01 8.550109e-02 0.0427505460
[53,] 9.706757e-01 5.864857e-02 0.0293242842
[54,] 9.605838e-01 7.883248e-02 0.0394162418
[55,] 9.637304e-01 7.253918e-02 0.0362695912
[56,] 9.590986e-01 8.180280e-02 0.0409014009
[57,] 9.427270e-01 1.145460e-01 0.0572730191
[58,] 9.415923e-01 1.168155e-01 0.0584077315
[59,] 9.857264e-01 2.854718e-02 0.0142735876
[60,] 9.860272e-01 2.794564e-02 0.0139728204
[61,] 9.870213e-01 2.595730e-02 0.0129786502
[62,] 9.924233e-01 1.515350e-02 0.0075767493
[63,] 9.818901e-01 3.621989e-02 0.0181099445
[64,] 9.757558e-01 4.848833e-02 0.0242441650
[65,] 9.946311e-01 1.073784e-02 0.0053689216
[66,] 9.968310e-01 6.337955e-03 0.0031689773
[67,] 9.992412e-01 1.517555e-03 0.0007587773
> postscript(file="/var/fisher/rcomp/tmp/1goh21355564548.ps",horizontal=F,onefile=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/fisher/rcomp/tmp/24x871355564548.ps",horizontal=F,onefile=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/fisher/rcomp/tmp/3058o1355564548.ps",horizontal=F,onefile=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/fisher/rcomp/tmp/4oxu41355564548.ps",horizontal=F,onefile=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/fisher/rcomp/tmp/5mksn1355564548.ps",horizontal=F,onefile=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 = 82
Frequency = 1
1 2 3 4 5 6
-3.9993428 -7.5146297 -8.9911896 -7.7616932 -3.9339203 -3.2066281
7 8 9 10 11 12
-0.8963015 -7.1004931 1.7680950 0.1860539 -4.5573397 -7.6295575
13 14 15 16 17 18
-6.6920061 0.1432400 -0.6989931 -0.2931115 2.1976392 4.7201319
19 20 21 22 23 24
2.6603599 -4.8273308 -6.9072257 -2.8163787 -12.3189262 -6.3298009
25 26 27 28 29 30
-4.3402201 -9.2744085 2.2645699 -5.0077526 -8.7736398 -3.1159937
31 32 33 34 35 36
-12.2409515 -8.1013669 -7.4954517 -1.4501269 -7.5278324 -7.4941640
37 38 39 40 41 42
-6.1358574 -11.6810532 -5.0558139 -5.4499425 -7.2638050 -2.1934530
43 44 45 46 47 48
7.4682342 12.1309316 11.4568846 2.0216582 7.9635901 3.5105151
49 50 51 52 53 54
2.6902866 5.7111080 3.5182830 7.8797823 3.2344275 2.0538863
55 56 57 58 59 60
10.3827552 11.8238849 8.9044783 4.6857854 12.7742954 13.4963756
61 62 63 64 65 66
7.2543586 5.7678266 9.1097734 7.7915181 11.3077982 11.5324971
67 68 69 70 71 72
6.5118823 -5.9716009 -3.3462418 -2.1060895 -8.0900656 9.0152247
73 74 75 76 77 78
-5.6595330 -4.1888340 6.4499060 5.2115457 7.6432809 4.5040865
79 80 81 82
1.4024938 -0.3548559 1.6717226 -2.0272448
> postscript(file="/var/fisher/rcomp/tmp/644mr1355564548.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> dum <- cbind(lag(myerror,k=1),myerror)
> dum
Time Series:
Start = 0
End = 82
Frequency = 1
lag(myerror, k = 1) myerror
0 -3.9993428 NA
1 -7.5146297 -3.9993428
2 -8.9911896 -7.5146297
3 -7.7616932 -8.9911896
4 -3.9339203 -7.7616932
5 -3.2066281 -3.9339203
6 -0.8963015 -3.2066281
7 -7.1004931 -0.8963015
8 1.7680950 -7.1004931
9 0.1860539 1.7680950
10 -4.5573397 0.1860539
11 -7.6295575 -4.5573397
12 -6.6920061 -7.6295575
13 0.1432400 -6.6920061
14 -0.6989931 0.1432400
15 -0.2931115 -0.6989931
16 2.1976392 -0.2931115
17 4.7201319 2.1976392
18 2.6603599 4.7201319
19 -4.8273308 2.6603599
20 -6.9072257 -4.8273308
21 -2.8163787 -6.9072257
22 -12.3189262 -2.8163787
23 -6.3298009 -12.3189262
24 -4.3402201 -6.3298009
25 -9.2744085 -4.3402201
26 2.2645699 -9.2744085
27 -5.0077526 2.2645699
28 -8.7736398 -5.0077526
29 -3.1159937 -8.7736398
30 -12.2409515 -3.1159937
31 -8.1013669 -12.2409515
32 -7.4954517 -8.1013669
33 -1.4501269 -7.4954517
34 -7.5278324 -1.4501269
35 -7.4941640 -7.5278324
36 -6.1358574 -7.4941640
37 -11.6810532 -6.1358574
38 -5.0558139 -11.6810532
39 -5.4499425 -5.0558139
40 -7.2638050 -5.4499425
41 -2.1934530 -7.2638050
42 7.4682342 -2.1934530
43 12.1309316 7.4682342
44 11.4568846 12.1309316
45 2.0216582 11.4568846
46 7.9635901 2.0216582
47 3.5105151 7.9635901
48 2.6902866 3.5105151
49 5.7111080 2.6902866
50 3.5182830 5.7111080
51 7.8797823 3.5182830
52 3.2344275 7.8797823
53 2.0538863 3.2344275
54 10.3827552 2.0538863
55 11.8238849 10.3827552
56 8.9044783 11.8238849
57 4.6857854 8.9044783
58 12.7742954 4.6857854
59 13.4963756 12.7742954
60 7.2543586 13.4963756
61 5.7678266 7.2543586
62 9.1097734 5.7678266
63 7.7915181 9.1097734
64 11.3077982 7.7915181
65 11.5324971 11.3077982
66 6.5118823 11.5324971
67 -5.9716009 6.5118823
68 -3.3462418 -5.9716009
69 -2.1060895 -3.3462418
70 -8.0900656 -2.1060895
71 9.0152247 -8.0900656
72 -5.6595330 9.0152247
73 -4.1888340 -5.6595330
74 6.4499060 -4.1888340
75 5.2115457 6.4499060
76 7.6432809 5.2115457
77 4.5040865 7.6432809
78 1.4024938 4.5040865
79 -0.3548559 1.4024938
80 1.6717226 -0.3548559
81 -2.0272448 1.6717226
82 NA -2.0272448
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -7.5146297 -3.9993428
[2,] -8.9911896 -7.5146297
[3,] -7.7616932 -8.9911896
[4,] -3.9339203 -7.7616932
[5,] -3.2066281 -3.9339203
[6,] -0.8963015 -3.2066281
[7,] -7.1004931 -0.8963015
[8,] 1.7680950 -7.1004931
[9,] 0.1860539 1.7680950
[10,] -4.5573397 0.1860539
[11,] -7.6295575 -4.5573397
[12,] -6.6920061 -7.6295575
[13,] 0.1432400 -6.6920061
[14,] -0.6989931 0.1432400
[15,] -0.2931115 -0.6989931
[16,] 2.1976392 -0.2931115
[17,] 4.7201319 2.1976392
[18,] 2.6603599 4.7201319
[19,] -4.8273308 2.6603599
[20,] -6.9072257 -4.8273308
[21,] -2.8163787 -6.9072257
[22,] -12.3189262 -2.8163787
[23,] -6.3298009 -12.3189262
[24,] -4.3402201 -6.3298009
[25,] -9.2744085 -4.3402201
[26,] 2.2645699 -9.2744085
[27,] -5.0077526 2.2645699
[28,] -8.7736398 -5.0077526
[29,] -3.1159937 -8.7736398
[30,] -12.2409515 -3.1159937
[31,] -8.1013669 -12.2409515
[32,] -7.4954517 -8.1013669
[33,] -1.4501269 -7.4954517
[34,] -7.5278324 -1.4501269
[35,] -7.4941640 -7.5278324
[36,] -6.1358574 -7.4941640
[37,] -11.6810532 -6.1358574
[38,] -5.0558139 -11.6810532
[39,] -5.4499425 -5.0558139
[40,] -7.2638050 -5.4499425
[41,] -2.1934530 -7.2638050
[42,] 7.4682342 -2.1934530
[43,] 12.1309316 7.4682342
[44,] 11.4568846 12.1309316
[45,] 2.0216582 11.4568846
[46,] 7.9635901 2.0216582
[47,] 3.5105151 7.9635901
[48,] 2.6902866 3.5105151
[49,] 5.7111080 2.6902866
[50,] 3.5182830 5.7111080
[51,] 7.8797823 3.5182830
[52,] 3.2344275 7.8797823
[53,] 2.0538863 3.2344275
[54,] 10.3827552 2.0538863
[55,] 11.8238849 10.3827552
[56,] 8.9044783 11.8238849
[57,] 4.6857854 8.9044783
[58,] 12.7742954 4.6857854
[59,] 13.4963756 12.7742954
[60,] 7.2543586 13.4963756
[61,] 5.7678266 7.2543586
[62,] 9.1097734 5.7678266
[63,] 7.7915181 9.1097734
[64,] 11.3077982 7.7915181
[65,] 11.5324971 11.3077982
[66,] 6.5118823 11.5324971
[67,] -5.9716009 6.5118823
[68,] -3.3462418 -5.9716009
[69,] -2.1060895 -3.3462418
[70,] -8.0900656 -2.1060895
[71,] 9.0152247 -8.0900656
[72,] -5.6595330 9.0152247
[73,] -4.1888340 -5.6595330
[74,] 6.4499060 -4.1888340
[75,] 5.2115457 6.4499060
[76,] 7.6432809 5.2115457
[77,] 4.5040865 7.6432809
[78,] 1.4024938 4.5040865
[79,] -0.3548559 1.4024938
[80,] 1.6717226 -0.3548559
[81,] -2.0272448 1.6717226
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -7.5146297 -3.9993428
2 -8.9911896 -7.5146297
3 -7.7616932 -8.9911896
4 -3.9339203 -7.7616932
5 -3.2066281 -3.9339203
6 -0.8963015 -3.2066281
7 -7.1004931 -0.8963015
8 1.7680950 -7.1004931
9 0.1860539 1.7680950
10 -4.5573397 0.1860539
11 -7.6295575 -4.5573397
12 -6.6920061 -7.6295575
13 0.1432400 -6.6920061
14 -0.6989931 0.1432400
15 -0.2931115 -0.6989931
16 2.1976392 -0.2931115
17 4.7201319 2.1976392
18 2.6603599 4.7201319
19 -4.8273308 2.6603599
20 -6.9072257 -4.8273308
21 -2.8163787 -6.9072257
22 -12.3189262 -2.8163787
23 -6.3298009 -12.3189262
24 -4.3402201 -6.3298009
25 -9.2744085 -4.3402201
26 2.2645699 -9.2744085
27 -5.0077526 2.2645699
28 -8.7736398 -5.0077526
29 -3.1159937 -8.7736398
30 -12.2409515 -3.1159937
31 -8.1013669 -12.2409515
32 -7.4954517 -8.1013669
33 -1.4501269 -7.4954517
34 -7.5278324 -1.4501269
35 -7.4941640 -7.5278324
36 -6.1358574 -7.4941640
37 -11.6810532 -6.1358574
38 -5.0558139 -11.6810532
39 -5.4499425 -5.0558139
40 -7.2638050 -5.4499425
41 -2.1934530 -7.2638050
42 7.4682342 -2.1934530
43 12.1309316 7.4682342
44 11.4568846 12.1309316
45 2.0216582 11.4568846
46 7.9635901 2.0216582
47 3.5105151 7.9635901
48 2.6902866 3.5105151
49 5.7111080 2.6902866
50 3.5182830 5.7111080
51 7.8797823 3.5182830
52 3.2344275 7.8797823
53 2.0538863 3.2344275
54 10.3827552 2.0538863
55 11.8238849 10.3827552
56 8.9044783 11.8238849
57 4.6857854 8.9044783
58 12.7742954 4.6857854
59 13.4963756 12.7742954
60 7.2543586 13.4963756
61 5.7678266 7.2543586
62 9.1097734 5.7678266
63 7.7915181 9.1097734
64 11.3077982 7.7915181
65 11.5324971 11.3077982
66 6.5118823 11.5324971
67 -5.9716009 6.5118823
68 -3.3462418 -5.9716009
69 -2.1060895 -3.3462418
70 -8.0900656 -2.1060895
71 9.0152247 -8.0900656
72 -5.6595330 9.0152247
73 -4.1888340 -5.6595330
74 6.4499060 -4.1888340
75 5.2115457 6.4499060
76 7.6432809 5.2115457
77 4.5040865 7.6432809
78 1.4024938 4.5040865
79 -0.3548559 1.4024938
80 1.6717226 -0.3548559
81 -2.0272448 1.6717226
> 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/fisher/rcomp/tmp/7tvkb1355564548.ps",horizontal=F,onefile=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/fisher/rcomp/tmp/8qghy1355564548.ps",horizontal=F,onefile=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/fisher/rcomp/tmp/9f56o1355564548.ps",horizontal=F,onefile=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/fisher/rcomp/tmp/10tqul1355564548.ps",horizontal=F,onefile=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/fisher/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/fisher/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/fisher/rcomp/tmp/11zoa91355564548.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/fisher/rcomp/tmp/12t1jp1355564548.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/fisher/rcomp/tmp/13u9ud1355564548.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/fisher/rcomp/tmp/14ev9r1355564548.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/fisher/rcomp/tmp/158vrr1355564548.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/fisher/rcomp/tmp/16xo251355564548.tab")
+ }
>
> try(system("convert tmp/1goh21355564548.ps tmp/1goh21355564548.png",intern=TRUE))
character(0)
> try(system("convert tmp/24x871355564548.ps tmp/24x871355564548.png",intern=TRUE))
character(0)
> try(system("convert tmp/3058o1355564548.ps tmp/3058o1355564548.png",intern=TRUE))
character(0)
> try(system("convert tmp/4oxu41355564548.ps tmp/4oxu41355564548.png",intern=TRUE))
character(0)
> try(system("convert tmp/5mksn1355564548.ps tmp/5mksn1355564548.png",intern=TRUE))
character(0)
> try(system("convert tmp/644mr1355564548.ps tmp/644mr1355564548.png",intern=TRUE))
character(0)
> try(system("convert tmp/7tvkb1355564548.ps tmp/7tvkb1355564548.png",intern=TRUE))
character(0)
> try(system("convert tmp/8qghy1355564548.ps tmp/8qghy1355564548.png",intern=TRUE))
character(0)
> try(system("convert tmp/9f56o1355564548.ps tmp/9f56o1355564548.png",intern=TRUE))
character(0)
> try(system("convert tmp/10tqul1355564548.ps tmp/10tqul1355564548.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
6.300 1.611 7.904