*********************************************************************
*** Set the CALENDAR and ALLOCATE range. CAL to start in November, 1959:
**********************************************************************

calendar 1959 1 12
allocate 1969:12

**********************************************************************
*** Read in the Data:
**********************************************************************

open data basics.wks
data(format=wks,org=columns) / rate m1 m2 ip ppi
* For BASICS.RAT, you would use "data(for=rats) / rate m1 m2 ip ppi"
* For BASICS.XLS, you would use "data(for=xls,org=columns) / rate m1 m2 ip ppi"
* For BASICS.PRN, you would use "data(for=prn,org=columns) / rate m1 m2 ip ppi"
* For BASICS.DAT, you would use "data(for=free,org=columns) / rate m1 m2 ip ppi"

**********************************************************************
*** Look at the data numerically:
**********************************************************************

table
table / ppi rate

statistics rate

print / m1 m2
print(picture='*.#') / m1 m2

print 1960:1 1960:12 m1

**********************************************************************
*** Data transformations and new series
**********************************************************************

* Need first difference of the produce price index.
* The following two instructions are equivalent:
set ppidiff = ppi - ppi{1}
diff ppi / ppidiff

* Now, need (M1(t) - M1(t-3)):
set m1diff = m1 - m1{3}
sho mem
print(picture='*.##') / ppi ppidiff m1 m1diff

* And quarter-to-quarter and annual growth rates for M2 and PPI:

*  M2(t) - M2(t-1) divided by M2(t-1):
set grm2 = (m2 - m2{1})/m2{1}

*  PPI(t) - PPI(t-1) divided by PPI(t-1):
set grppi = (ppi - ppi{1})/ppi{1}

*  IP(t) - IP(t-1) divided by PPI(t-1):
set grip = 100*(ip - ip{1})/ip{1}
* Annualized growth rates:
set anngrppi = 100*( (ppi/ppi{1})**12 -1.0)
set anngrip  = 100*( (ip/ip{1})**12 -1.0)

* Need a three period weighted moving average:
set pratio = ppidiff/ppi
set ppisum = pratio + pratio{1} + pratio{2}

* Or:
set ppisum = (ppidiff/ppi) + (ppidiff{1}/ppi{1}) + (ppidiff{2}/ppi{2})

* Or:
set pratio = ppidiff/ppi
filter pratio / ppisum
# 1 2
# 1.0 1.0

* Simple trend and squared-trend series, commonly used:
*set trend = t
*set trendsq = t**2
**********************************************************************
*** Graph the data:
**********************************************************************

* Look at the rate and index series together:
graph(key=upleft) 3
# rate
# ip
# ppi


* Graph the money supply series
graph(key=below) 2
# m1
# m2

* Again, now using some labelling options:
graph(key=upleft,header='Real Money Supply (M1 and M2)',subhead='Billions US$, Seasonally Adjusted') 2
# m1
# m2

* Use SMPL to change default range:
smpl 1959:11 1960:12
graph(key=upleft) 3
# rate
# ip
# ppi

* Reset the SMPL back to full CALENDAR-ALLOCATE range:
smpl

* Some SCATTER plots:
scatter 1
# grppi rate


scatter(vlabel='Interest Rate',hlabel='PPI Growth Rate',header='Interest Rates vs. PPI Growth') 1
# grppi rate


**********************************************************************
*** Fit some least squares regressions.
**********************************************************************
* Example 4.2 from 3rd (1991) Edition:
* Data are slightly different, so results do not match text exactly
linreg rate 1960:2 *
# constant ip m1diff ppisum

* Example 4.2 from 4th (1998) edition:
linreg rate /
# constant ip grm2 grppi{1}

**********************************************************************
*** Hypothesis testing:
**********************************************************************

linreg rate
# constant rate{1 to 6}

exclude
# rate{4 to 6}

* Same thing using TEST:
test
# 5   6   7
# 0.0 0.0 0.0

* Test that RATE{1}==1.0 and RATE{2}==0.0
test
# 2 3
# 1.0 0.0


* Test whether coefficient on GRM2 is equal to coeff. on GRPPI{1}:
* Because RESTRICT test linear combinations, to test
* H0: beta1 == beta2
* We write the test as:
* H0: beta1 - beta2 == 0.0

linreg rate
# constant ip grm2 grppi{1}

* Test 1 resriction: 1.0*(coefficient 3) - 1.0*(coefficient 4) == 0.0
restrict 1
# 3     4
# 1.0  -1.0  0.0


**********************************************************************
*** Forecasting:
**********************************************************************
* Using OLS model:

linreg rate 1960:1 *
# constant ip grm2 grppi{1}

* Fitted values:
prj fitted

graph(key=below,header='Actual vs. Fitted') 2
# rate  1960:1 *
# fitted 1960:1 *


* Forecasts (same results as PRJ for static models like this one)

linreg(define=irateeq) rate 1960:1 *
# constant ip grm2 grppi{1}

forecast 1 14 1968:10
# irateeq olsfore

graph(header='Interest Rate and Forecast',key=below) 2
# rate
# olsfore

**********************************************************************
*** Scalar and matrix examples:
**********************************************************************

compute a = 1.0
display a

make xmat 1960:1 *
# constant ip grm2 grppi{1}
make ymat 1960:1 *
# rate

compute invXX = inv(tr(xmat)*(xmat))
compute XY    = tr(xmat)*ymat
compute beta = invXX*XY
display beta

**********************************************************************
*** Estimate with Cochrane-Orcutt type serial correlation correction:
**********************************************************************

* Example 6.5 in 3rd edition:
ar1(method=corc) rate 1960:1 *
# constant ip m1 ppisum

* Example 6.6 in 4th edition:
smpl 1960:1 *
linreg rate / olsres
# constant ip grm2 grppi{1}

ar1(method=corc) rate / ar1res
# constant ip grm2 grppi{1}

graph(key=upleft) 2
# olsres
# ar1res

smpl


* Forecast using AR1 model (Example 8.3 in in P&R 4th Edition):
ar1(method=corc,define=ar1eq) rate 1960:1 *
# constant ip grm2 grppi{1}

smpl 1969:1 1969:12
forecast 1
# ar1eq ar1fore
* Compare with actual and "OLSFORE" series computed from OLS
*  model above:
graph(header='Interest Rate and Forecasts',key=below) 3
# rate
# ar1fore
# olsfore

smpl

**********************************************************************
*** Display VCV Matrix:
**********************************************************************

linreg(vcv) rate 1960:1 *
# constant ip grm2 grppi{1}


**********************************************************************
*** Estimate using NLLS:
**********************************************************************

nonlin b0 b1 b2 b3
frml f1 rate = b0 + b1*ip + b2*m1diff + b3*ppisum
nlls(frml=f1) rate 1960:2 1980:12

frml f2 rate = b0 + b1*ip + b2*grm2 + b3*grppi{1}
nlls(frml=f2) rate 1960:1 1995:8