Using gret l for Principles of Econometrics, 4th Edition
FGLS and Other Estimators
The feasible GLS estimator of the AR(p) model can be estimated using gretl in a number of ways. For first order autocorrelated models the ar1 command can be used. There are a num
ber of estimators available by option including the Cochrane-Orcutt (iterated), the Prais-Winsten (iterated), and the Hildreth-Lu search procedure. Examples are:
1 list x = d_u const
2 ar1 inf x # Cochrane-Orcutt (default)
3 ar1 inf x —pwe # Prais-Winsten
4 ar1 inf x —hilu —no-corc # Hildreth-Lu
The results are collected in a model table below.
AR(1) Errors
Dependent variable: inf
(CO) |
(PW) |
(HL) |
|
const |
0.7609** |
0.7862** |
0.7608** |
(0.1238) |
(0.1218) |
(0.1245) |
|
d_u |
-0.6944** |
-0.7024** |
-0.6953* |
(0.2429) |
(0.2430) |
(0.2430) |
|
P |
0.55739 |
0.55825 |
.56 |
n |
89 |
90 |
89 |
R2 |
0.3407 |
0.3418 |
0.3406 |
Standard errors in parentheses
* indicates significance at the 10 percent level
** indicates significance at the 5 percent level
CO = Cochrane Orcutt, PW=Prais-Winsten, HL=Hildreth-Lu
You can see that there are minor differences produced by these options. If the —no-corc option is not used with --hilu then the Hildreth-Lu estimator is modified slightly to perform additional iterations as the end. Notice that the Prais-Winsten is the only procedure to use all 90 observations.
For higher order models there are two commands worth taking note of. The ar command estimates a linear regression with arbitrary autocorrelation structure. It uses a generalization of the Cochrane-Orcutt iterative procedure to obtain estimates.
The other estimator is arima, the syntax for which appears below:
p d q [ ; P D Q ] ; depvar [ indepvars ]
—verbose [print details of iterations)
—vcv (print covariance matrix)
--nessian (see below)
—opg (see below)
—nc (do not include a constant)
—conditional (use conditional maximum likelihood) —x-12-arima (use X-12-ARIMA for estimation)
— lbf gs (use L-BFGS-B maximizer)
—y-dif f-only (ARIMAX special, see below) --save-enat (see below) arima 1 0 2 ; у
arima 202 ; у 0 xl x2 --verbose arima Oil; Oil; у —nc
The default estimation method for arima in gretl is to estimate the parameters of the model using the “native” gretl ARMA functionality, with estimation by exact maximum likelihood using the Kalman filter. You can estimate the parameters via conditional maximum likelihood as well.
Estimating the simple AR(1) regression using these estimators is done:
1 ar 1 ; inf x
2 arima 100; inf x
For the ar command, list the lag numbers for the desired residuals. In the case of AR(1) this is just 1. This is followed by a semicolon and then the regression to estimate. The arima syntax is similar, except you specify p, d, and q, where p is the order of the desired autocorrelation, d is the number of differences to take of the time-series, and q is the order of any moving average terms you might have in the residuals.
The outcome for the simple ARIMA(1,0,0) ia
ARMAX, using observations 1987:2-2009:3 (T = 90)
Estimated using Kalman filter (exact ML)
Dependent variable: inf Standard errors based on Hessian
coefficient std. error z p-value
|
S. D. dependent var |
0.636819 |
S. D. of innovations |
0.510937 |
Akaike criterion |
142.9118 |
Hannan-Quinn |
146.9441 |
maginary Modulus |
Frequency |
0.0000 1.7895 |
0.0000 |
Mean dependent var 0.791111 Mean of innovations -0.003996 Log-likelihood -67.45590 Schwarz criterion 152.9110 Real |
AR Root 1 1.7895 |
These are very similar to the ones above. The coefficient labeled phi_1 is the estimate of the autocorrelation parameter. The root of this equation is 1/phi_1. The roots (or modulus) must be greater than 1 in absolute value in order for the model to be stationary.
1 open "@gretldirdatapoeokun. gdt"
2 set echo off
3 # change variable attributes
4 setinfo g - d "percentage change in U. S. Gross Domestic Product, seasonally
5 adjusted" -n "Real GDP growth"
6 setinfo u - d "U. S. Civilian Unemployment Rate (Seasonally adjusted)" - n
7 "Unemployment Rate"
8
8 # plot series and save output to files
9 gnuplot g —with-lines —time-series —output="@workdirokun_g. plt"
10 gnuplot u —with-lines —time-series —output="@workdirokun_u. plt"
12
11 # graphing multiple time-series
12 scatters g u --with-lines
15
13 diff u
14 setinfo d_u - d "Change in U. S. Civilian Unemployment
15 Rate (Seasonally adjusted)" - n
16 "D. Unemployment Rate"
17 scatters g d_u --with-lines --output=display
21
18 # distributed lag models
19 ols d_u const g(0 to -3)
20 smpl 1986:1 2009:3
21 ols d_u const g(0 to -2)
26
22 gnuplot g g_1
28
23 # correlogram and confidence interval
24 corrgm g 12
31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 |
matrix ac = corrgm(g, 12) matrix lb = ac[,1]-1.96/sqrt($nobs) matrix ub = ac[,1]+1.96/sqrt($nobs) matrix all = lb~ac[,1]~ub colnames(all, "Lower AC Upper ")
printf "nAutocorrelations and 95%% confidence intervalsn %9.4fn", all
# Phillips curve
open "@gretldirdatapoephillips_aus. gdt" diff u
setinfo inf - d "Australian Inflation Rate" - n "Inflation Rate" setinfo d_u - d "Change in Australian Civilian
Unemployment Rate (Seasonally adjusted)" - n
"D. Unemployment Rate" scatters inf d_u --with-lines
ols inf const d_u series ehat = $uhat gnuplot ehat —time-series corrgm ehat
# LM tests
ols ehat const d_u ehat(-1) scalar NR2 = $trsq pvalue X 1 NR2
ols ehat const d_u ehat(-1 to -4) scalar NR2 = $trsq pvalue X 4 NR2
ols inf const d_u modtest 1 —autocorr modtest 4 —autocorr —quiet
# HAC standard errors
open "@gretldirdatapoephillips_aus. gdt" set hac_kernel bartlett set hac_lag nw2 diff u
ols inf const d_u modeltab add
ols inf const d_u --robust modeltab add modeltab show modeltab free
# nonlinear least squares estimation of regression w/AR(1) errors open "@gretldirdatapoephillips_aus. gdt"
diff u
ols inf const d_u --quiet
82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 |
scalar betal = $coeff(const) scalar beta2 = $coeff(d_u) scalar rho = 0
nls inf = beta1*(1-rho) + rho*inf(-1) + beta2*(d_u-rho*d_u(-1)) params rho beta1 beta2 end nls
scalar delta = $coeff(beta1)*(1-$coeff(rho)) scalar delta1 = -$coeff(rho)*$coeff(beta2)
printf "nThe estimated delta is %.3f and the estimated delta1 is %.3f.n",delta, delta1 scalar sser=$ess
# estimation of more general model ols inf const inf(-1) d_u(0 to -1) scalar sseu=$ess
scalar fstat = (sser-sseu)/(sseu/$df) pvalue X 1 fstat pvalue F 1 $df fstat omit d_u(-1)
ols inf const inf(-1) d_u(0 to -1) modeltab add
ols inf const inf(-1) d_u(0) modeltab add modeltab show modeltab free
# model selection function
function matrix modelsel (series y, list xvars) ols y xvars —quiet scalar sse = $ess scalar N = $nobs scalar K = nelem(xvars) scalar aic = ln(sse/N)+2*K/N scalar bic = ln(sse/N)+K*ln(N)/N matrix A = { K, N, aic, bic} printf "nRegressors: %sn",varname(xvars)
printf "K = %d, N = %d, AIC = %.4f SC = %.4f.n",K, N,aic, bic return A end function
# using the modelsel function
list x = const inf(-1) d_u(0 to -1)
matrix a = modelsel(inf, x)
list x0 = const
matrix b = modelsel(inf, x)
list x = const d_u inf(-1)
# putting the model selection results into a matrix open "@gretldirdatapoephillips_aus. gdt"
133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 |
diff u
smpl 1988:3 2009:3 matrix A = {} scalar q = 0 loop p = 1..6 —quiet if p = 1
list x = const inf(-1) d_u
else
list x = const inf(-1 to -$p) d_u endif
matrix a = $p~q~modelsel(inf, x) matrix A = A | a modelsel(inf, x) endloop scalar q = 1 loop p = 1..6 —quiet if p = 1
list x = const inf(-1) d_u(0 to -1)
else
list x = const inf(-1 to -$p) d_u(0 to -1) endif
matrix a = $p~q~modelsel(inf, x) matrix A = A | a endloop
colnames(A,"p q K N AIC SC ") print A
smpl full
ols inf const inf(-1 to -4) d_u —robust
# improved modelsel2 function for ARDL function matrix modelsel2 (list xvars)
ols xvars --quiet
scalar sse = $ess
scalar N = $nobs
scalar K = nelem(xvars)-1
scalar aic = ln(sse/N)+2*K/N
scalar bic = ln(sse/N)+K*ln(N)/N
matrix A = { K, N, aic, bic}
# printf "nDependent variable and Regressors: %sn",varname(xvars)
# printf "K = %d, N = %d, AIC = %.4f SC = %.4f.n",K, N,aic, bic return A
end function
# using modelsel2
open "@gretldirdatapoeokun. gdt" diff u
smpl 1986:1 2009:3
matrix A = {}
loop p = 0..2 —quiet
loop q = 1..3 —quiet
184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 |
if p=0
list vars = d_u g(0 to - q) const
else
list vars = d_u(0 to - p) g(0 to - q) const endif
matrix a = p~q~modelsel2(vars) matrix A = A | a endloop endloop
colnames(A,"p q K N AIC SC ") print A
function modelsel clear smpl full
ols d_u(0 to -1) g(0 to -1) const loop i=1..4
modtest $i —autocorr —quiet endloop
open "@gretldirdatapoeokun. gdt" smpl 1986:3 2009:3 matrix A = {} scalar q=0
loop p = 1..5 —quiet
list vars = g(0 to - p) const matrix a = p~q~modelsel2(vars) matrix A = A | a endloop
colnames(A,"p q K N AIC SC ") print A
function modelsel clear
# loop to test for autocorrelation in ARDL open "@gretldirdatapoephillips_aus. gdt" diff u
ols inf(0 to -1) d_u const loop i=1..5
modtest $i —autocorr —quiet endloop
# loop to test for autocorrelation at several lags open "@gretldirdatapoeokun. gdt"
ols g(0 to -2) const series res = $uhat corrgm res loop i = 1..4
modtest $i —autocorr —quiet endloop
# model selection for Okun data open "@gretldirdatapoeokun. gdt"
235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 |
smpl 1986:3 2009:3 matrix A = {} scalar q=0
loop p = 1..5 —quiet
list vars = g(0 to - p) const matrix a = p~q~modelsel2(vars) matrix A = A | a endloop
colnames(A,"p q K N AIC SC ") print A
# estimation of preferred model and a forecast open "@gretldirdatapoeokun. gdt"
ols g(0 to -2) const dataset addobs 3
fcast 2009:4 2010:2 —plot="@workdirar2plot1.plt"
# multiplier analysis
open "@gretldirdatapoeokun. gdt"
matrix y = { g }
scalar T = $nobs
matrix sm1 = zeros(T,1)
scalar a = .38
smpl 1 round((T+1)/2)
scalar stv = mean(y)
smpl full
loop i=1..T --quiet if i = 1
matrix sm1[i]=stv
else
matrix sm1[$i]=a*y[$i] + (1-a)*sm1[i-1] endif endloop
series exsm = sm1
gnuplot g exsm —time-series
scalar a = .8 loop i=1..T --quiet if i = 1
matrix sm1[i]=stv
else
matrix sm1[$i]=a*y[$i] + (1-a)*sm1[i-1] endif endloop
series exsm8 = sm1
gnuplot g exsm8 —time-series
open "@gretldirdatapoeokun. gdt" diff u
ols d_u(0 to -1) g(0 to -1) const scalar b0 = $coeff(g)
286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 |
scalar b1 =$coeff(d_u_1)*b0+$coeff(g_1) scalar b2 = b1*$coeff(d_u_1) scalar b3 = b2*$coeff(d_u_1)
# Matrix & Series Plot
open "@gretldirdatapoeokun. gdt" diff u
ols d_u(0 to -1) g(0 to -1) const scalar h = 8
matrix mult = zeros(h,2) loop i=1..h
mult[i,1] = i-1 scalar b0 = $coeff(g)
scalar b1 = $coeff(d_u_1)*b0+$coeff(g_1) if i=1
mult[i,2]=b0 elif i=2
mult[i,2]=b1
else
mult[i,2]=mult[i-1,2]*$coeff(d_u_1)
endif
endloop
gnuplot 2 1 —matrix=mult —output=display —with-lines —suppress-fitted
printf "nThe impact and delay multipliers are n %10.5fn", mult
nulldata 8 --preserve series m = mult[,2] series lag = mult[,1]
setinfo m - d "Multipliers" - n "Multiplier"
gnuplot m index —with-lines —output=display —suppress-fitted
# appendix
open "@gretldirdatapoephillips_aus. gdt" diff u
setinfo inf - d "Australian Inflation Rate" - n "Inflation Rate" setinfo d_u - d "Change in Australian Civilian
Unemployment Rate (Seasonally adjusted)" - n
"D. Unemployment Rate"
# Durbin-Watson with p-value list x = d_u const
ols inf x
scalar dw_p = $dwpval print dw_p
# various ways to estimate AR(1) regression ar1 inf x
modeltab add ar1 inf x --pwe
337 modeltab add
338 ar1 inf x —hilu —no-corc
339 modeltab add
340 modeltab show
341 modeltab free
342
342 ar 1 ; inf x
343 arima 10 0; inf x
residual |
0 2 4 6 8 to 12 lag |
Breusch-Godfrey teat for autocorrelation up to order 4 OLS, using observations 1987:2-2009:3 (T = 90) Dependent variable: uhat
COn3t
d_u uhat_l uhat_2 uhat_3 l. lll 4
Unadjusted R-squared = 0.407466
Test statistic: LMF - 14.440976,
with P-value = P(F(4,84) > 14.441) = S. lSe-009
Alternative statistic: TR~2 = 36.671897,
with p-value = P(Chi-square(4) > 36.6719) = 2.ІЄ-007
Ljung-Box Q1 = 82.4327,
with p-value = P(Chi-square(4) > 82.4327) = 5.3ІЄ-017
Using numerical derivatives
Tolerance = 1.81899e-012
Convergence achieved after 21 iterations
Model 2: NL5, using observations 1987:3-2009:3 (T = 89) inf = betal*(1-rho) + rho*inf(-l) + beta2*<d_u-rho*d_u(-1))
estimate std. error t-ratio p-value
Figure 9.13: Nonlinear least squares results for the AR(1) regression model.