Matlab

An Introduction to Matlab for Econometrics

John C. Frain

TEP Working Paper No. 0110 February 2010

Trinity Economics Papers
Department of Economics Trinity College Dublin

An Introduction to MATLAB for Econometrics
John C. Frain. February 2010
∗

Abstract This paper is an introduction to MATLAB for econometrics. It describes the MATLAB Desktop, contains a sample MATLAB session showing elementary MATLAB operations, gives details of data input/output, decision and loop structures, elementary plots, describes the LeSage econometrics toolbox and maximum likelihood using the LeSage toolbox. Various worked examples of the use of MATLAB in econometrics are also given. After reading this document the reader should be able to make better use of the MATLAB on-line help and manuals.

Contents
1 Introduction 1.1 1.2 Preliminaries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . The MATLAB Desktop . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.1 1.2.2 1.2.3 1.2.4 1.2.5 1.2.6 1.2.7 1.2.8 1.2.9
∗ Comments

4 4 6 6 7 8 8 9 9 9

The Command Window . . . . . . . . . . . . . . . . . . . . . . . . The Command History Window . . . . . . . . . . . . . . . . . . . The Start Button . . . . . . . . . . . . . . . . . . . . . . . . . . . . The Edit Debug window . . . . . . . . . . . . . . . . . . . . . . . . The Figure Windows . . . . . . . . . . . . . . . . . . . . . . . . . . The Workspace Browser . . . . . . . . . . . . . . . . . . . . . . . . The Help Browser . . . . . . . . . . . . . . . . . . . . . . . . . . .

The Path Browser . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 Miscellaneous Commands . . . . . . . . . . . . . . . . . . . . . . . 11

are welcome. My email address is frainj at tcd.ie

1

2 Vectors, Matrices and Arrays 2.1

11

A Sample MATLAB session . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.1.1 2.1.2 2.1.3 2.1.4 2.1.5 2.1.6 2.1.7 2.1.8 2.1.9 Entering Matrices . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 Basic Matrix operations . . . . . . . . . . . . . . . . . . . . . . . . 12 Kronecker Product . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 Examples of number formats . . . . . . . . . . . . . . . . . . . . . 15 fprintf function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 element by element operations . . . . . . . . . . . . . . . . . . . . 16 miscellaneous functions . . . . . . . . . . . . . . . . . . . . . . . . 17 sequences . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 Creating Special Matrices . . . . . . . . . . . . . . . . . . . . . . . 22

2.1.10 Random number generators . . . . . . . . . . . . . . . . . . . . . . 23 2.1.11 Extracting parts of a matrix, Joining matrices together to get a new larger matrix . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 2.1.12 Using sub-matrices on left hand side of assignment . . . . . . . . . 25 2.1.13 Stacking Matrices . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 2.1.14 Special Values . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 2.2 2.3 2.4 2.5 Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 Regression Example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 Simulation – Sample Size and OLS Estimates . . . . . . . . . . . . . . . . 31 Example – Macroeconomic Simulation with Matlab . . . . . . . . . . . . . 34 39

3 Data input/output 3.1 3.2 3.3 3.4 3.5 3.6 3.7

Native MatLab data ﬁles . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 Importing from Excel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 Reading from text ﬁles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 Exporting data to EXCEL, STATA and other programs . . . . . . . . . . 41 Stat/Transfer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 Formatted Output . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 Producing material for inclusion in a paper . . . . . . . . . . . . . . . . . 42 43 46

4 Decision and Loop Structures. 5 Elementary Plots

2

6 Systems of Regresssion Equations 6.1 6.2

47

Using Matlab to estimate systems of regression equations . . . . . . . . . 47 Exercise – Using Matlab to estimate a simultaneous equation systems . . 59 59 61 82 86

7 User written functions in MATLAB 8 The LeSage Econometric Toolbox 9 Maximum Likelihood Estimation using Numerical Techniques 10 Octave, Scilab and R

10.1 Octave . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 10.2 Scilab . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 10.3 R . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 A Functions etc. in LeSage Econometrics Toolbox 89

A.1 Regression . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 A.1.1 Programs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 A.1.2 Demonstrations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 A.1.3 Support functions . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 A.2 Utilities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 A.2.1 Utility Function Library . . . . . . . . . . . . . . . . . . . . . . . . 92 A.2.2 demonstration programs . . . . . . . . . . . . . . . . . . . . . . . . 94 A.3 Graphing Function Library . . . . . . . . . . . . . . . . . . . . . . . . . . 94 A.3.1 graphing programs . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 A.3.2 Demonstration Programs . . . . . . . . . . . . . . . . . . . . . . . 94 A.3.3 Support Functions . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 A.4 Regression Diagnostics Library . . . . . . . . . . . . . . . . . . . . . . . . 95 A.4.1 regression diagnostic programs . . . . . . . . . . . . . . . . . . . . 95 A.4.2 Demonstration Programs . . . . . . . . . . . . . . . . . . . . . . . 95 A.4.3 support functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 A.5 vector autoregressive function library . . . . . . . . . . . . . . . . . . . . . 96 A.5.1 VAR/BVAR functions . . . . . . . . . . . . . . . . . . . . . . . . . 96 A.5.2 Demonstration Programs . . . . . . . . . . . . . . . . . . . . . . . 96 A.5.3 Demonstration Programs - continued . . . . . . . . . . . . . . . . . 97 A.5.4 Support Functions . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 3

A.6 Co-integration Library . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 A.6.1 Co-integration testing routines . . . . . . . . . . . . . . . . . . . . 98 A.6.2 Demonstration Programs . . . . . . . . . . . . . . . . . . . . . . . 98 A.6.3 Support functions . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 A.7 Gibbs sampling convergence diagnostics functions . . . . . . . . . . . . . . 99 A.7.1 Convergence testing functions . . . . . . . . . . . . . . . . . . . . . 99 A.7.2 Demonstration Programs . . . . . . . . . . . . . . . . . . . . . . . 99 A.7.3 Support Functions . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 A.8 Distribution functions library . . . . . . . . . . . . . . . . . . . . . . . . . 100 A.8.1 pdf, cdf, inverse functions . . . . . . . . . . . . . . . . . . . . . . . 100 A.8.2 Random Samples . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 A.8.3 Demonstration and Test programs . . . . . . . . . . . . . . . . . . 102 A.8.4 Support Functions . . . . . . . . . . . . . . . . . . . . . . . . . . . 102 A.9 Optimization functions library . . . . . . . . . . . . . . . . . . . . . . . . 102 A.9.1 Optimization Functions . . . . . . . . . . . . . . . . . . . . . . . . 102 A.9.2 Demonstration Programs . . . . . . . . . . . . . . . . . . . . . . . 103 A.9.3 Support Functions . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 A.10 Spatial Econometrics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 A.10.1 Functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 A.10.2 Demonstration Programs . . . . . . . . . . . . . . . . . . . . . . . 104 A.10.3 Support Functions . . . . . . . . . . . . . . . . . . . . . . . . . . . 105

1
1.1

Introduction
Preliminaries

These notes are a guide for students of econometrics who wish to learn MATLAB in MS Windows. I deﬁne the functions of MATLAB using simple examples. To get the best beneﬁt from these notes you should read them sitting in front of a computer entering the various MATLAB instructions as you read the notes. The material in the ﬁrst three sections is elementary and will be required by all economists starting with MATLAB. The remaining sections contain some more advanced material and should be read as required. In these notes I have used a mono-spaced font for MATLAB instructions and for computer input and output. Descriptive material, explanations and commentary on the computer input/output is given in the current font. 4

While the ﬁrst aim of these notes is to get the reader started in the use of MATLAB for econometrics it should be pointed out that MATLAB has many uses in economics. In recent years it has been used widely in what has become known as computational economics. This area has applications in macroeconomics, determination of optimal policies and in ﬁnance. Recent references include Kendrick et al. (2006), Marimon and Scott (1999), Miranda and Fackler (2002) and Ljungqvist and Sargent (2004). I do not know of any book on MATLAB written speciﬁcally for economics. Creel (2008) is a set of lecture notes on econometrics which can be downloaded from the web. This contains examples of quantitative econometric analysis using GNU Octave which has a syntax similar to Matlab (see section 10.1). LeSage (1999) is a free econometrics toolbox available for download from http://www.spatial-econometrics.com/. This site also contains links to several other MATLAB resources useful in econometrics. A free ARCH/GARCH toolbox is available at http://http://www.kevinsheppard.com/ wiki/MFE_Toolbox. MathWorks, the composers of MATLAB have a list of books using MATLAB for Economics/Finance (see www.mathworks.com). They have also issued a new econometrics toolbox (see http://www.mathworks.com/products/econometrics/). The MathWorks overview of this toolbox indicates that is is targeted at econometric time series in ﬁnance. For advanced applications in applied probability Paolella (2006, 2007) are comprehensive accounts of computational aspects of probability theory using MATLAB. Higham and Higham (2005) is a good book on MATLAB intended for all users of MATLAB. Pratap (2006) is a good general “getting started” book. There are also many excellent books covering MATLAB for Engineers and/or Scientists which you might ﬁnd useful if you need to use MATLAB in greater depth. These notes can not give a comprehensive account of MATLAB. Your copy of MATLAB comes with one of the best on-line help systems available. Full versions of the manuals are available in portable document format on the web at http:/www.mathworks.com. MATLAB started life, in the late 70’s, as a computer program for handling matrix operations. Over the years it has been extended and the basic version of MATLAB now contains more than 1000 functions. Various “toolboxes” have also been written to add specialist functions to MATLAB. Anyone can extend MATLAB by adding their own functions and/or toolboxes. Any glance at an econometrics textbook shows that econometrics involves much matrix manipulation and MATLAB provides an excellent platform for implementing the various textbook procedures and other state of the art estimators. Before you use MATLAB to implement procedures from your textbook you must understand the matrix manipulations that are involved in the procedure. When you implement them you will understand the procedure better. Using a black box package may, in some cases, be easier but how often do you know exactly what the black box is producing. Using MATLAB for econometrics may appear to involve a lot of extra work but many students have found that it helps their understanding of both matrix theory and econometrics. In MATLAB as it all other packages it makes life much easier if you organise your work

5

properly. The procedure That I use is some variation of the following – 1. Set up a new directory for each project (e.g. s:\Matlab\project1) 2. Set up a shortcut for each project. The shortcut should specify that the program start in the data directory for the project. If all your work is on the same PC the shortcut is best stored on the desktop. If you are working on a PC in a computer lab you will not be able to use the desktop properly and the shortcut may be stored in the directory that you have set up for the project. If you have several projects in hand you should set up separate shortcuts and directories for each of them. Each shortcut should be renamed so that you can associate it with the relevant project. 3. Before starting MATLAB you are strongly advised to amend the options in Windows explorer so that full ﬁlenames (including any ﬁle extensions allocated to programs) appear in Windows Explorer and any other Windows ﬁle access menus.

1.2

The MATLAB Desktop

The MATLAB desktop has the following parts 1. The Command Window 2. The Command History Window 3. The Start Button 4. The Documents Window (including the Editor/(Debugger) and Array Editor 5. The Figure Windows 6. The Workspace Browser 7. The Help Browser 8. The Path Browser 1.2.1 The Command Window

The simplest use of the command window is as a calculator. With a little practice it may be as easy, if not easier, to use than a spreadsheet. Most calculations are entered almost exactly as one would write them.

>> 2+2 ans = 4 >> 3*2 ans = 6 6

The object ans contains the result of the last calculation of this kind. You may also create an object a which can hold the result of your calculation. >> a=3^3 a = 27 >> a a = 27 >> b=4^2+1 b = 17 >> b=4^2+1; % continuation lines >> 3+3 ... +3 ans = 9 Type each instruction in the command window, press enter and watch the answer. Note • The arithmetic symbols +, -, *, / and ^ have their usual meanings • The assignment operator = • the MATLAB command prompt >> • A ; at the end of a command does not produce output but the assignment is made or the command is completed • If a statement will not ﬁt on one line and you wish to continue it to a second type an ellipsis (. . . ) at the end of the line to be continued.

Individual instructions can be gathered together in an m-ﬁle and may be run together from that ﬁle (or script). An example of a simple m-ﬁle is given in the description of the Edit Debug window below. You may extend MATLAB by composing new MATLAB instructions using existing instructions gathered together in a script ﬁle. You may use the up down arrow keys to recall previous commands (from the current or earlier sessions) to the Command Window. You may the edit the recalled command before running it. Further access to previous commands is available through the command window. 1.2.2 The Command History Window

If you now look at the Command History Window you will see that as each command was entered it was copied to the Command History Window. This contains all commands 7

previously issued unless they are speciﬁcally deleted. To execute any command in the command history double click it with the left mouse button. To delete a commands from the history select them, right click the selection and select delete from the drop down menu. 1.2.3 The Start Button

This allows one to access various MATLAB functions and works in a manner similar to the Start button in Windows 1.2.4 The Edit Debug window

Clearly MATLAB would not be of much use if one was required, every time you used MATLAB, to enter your commands one by one in the Command Window. You can save your commands in an m-ﬁle and run the entire set of commands in the ﬁle. MATLAB also has facilities for changing the commands in the ﬁle, for deleting command or adding new commands to the ﬁle before running them. Set up and run the simple example below. We shall be using more elaborate examples later The Edit Window may be used to setup and edit M-ﬁles. Use File|New|m-file to open a new m-ﬁle. Enter the following in the ﬁle \vol\_sphere.m % vol_sphere.m % John C Frain revised 12 November 2006 % This is a comment line % This M-file calculates the volume of a sphere echo off r=2 volume = (4/3) * pi * r^3; string=[’The volume of a sphere of radius ’ ... num2str(r) ’ is ’ num2str(volume)]; disp(string) % change the value of r and run again Now Use File|Save As vol sphere.m. (This will be saved in your default directory if you have set up things properly check that this is working properly). Now return to the Command Window and enter vol sphere. If you have followed the instructions properly MATLAB will process this as if it were a MATLAB instruction. The edit window is a programming text editor with various features colour coded. Comments are in green, variables and numbers in black, incomplete character strings in red and language key-words in blue. This colour coding helps to identify errors in a program. The Edit window also provides debug features for use in ﬁnding errors and verifying programs. Additional features are available in the help ﬁles. 8

1.2.5

The Figure Windows

This is used to display graphics generated in MATLAB. Details will be given later when we are dealing with graphics. 1.2.6 The Workspace Browser

This is an option in the upper left hand window. Ensure that it is open by selecting the workspace tab. Compare this with the material in the command window. Note that it contains a list of the variables already deﬁned. Double clicking on an item in the workspace browser allows one to give it a new value. The contents of the workspace can also be listed by the whos command 1.2.7 The Help Browser

The Help Browser can be started by selecting the [?] icon from the desktop toolbar or by typing helpdesk or helpwin in the Command Window. You should look at the Overview and the Getting Help entries. There are also several demos which provide answers to many questions. MATLAB also has various built-in demos. To run these type demo at the command prompt or select demos from the start button The on-line documentation for MatLab is very good. The MatLab www site (http: /www.mathworks.com gives access to various MatLab manuals in pdf format. These may be downloaded and printed if required. (Note that some of these documents are very large and in many cases the on-line help is more accessible and is clearer. One can also type help at a command prompt to get a list of help topics in the Command Window. Then type help topic and details will be displayed in the command window. If, for example, you wish to ﬁnd a help ﬁle for the inverse of a matrix the command help inverse will not help as there is no function called inverse. In such a case one may enter lookfor inverse and the response will be

INVHILB Inverse Hilbert matrix. IPERMUTE Inverse permute array dimensions. ACOS ACOSD ACOSH ACOT ACOTD ACOTH ACSC ACSCD Inverse cosine, result in radians. Inverse cosine, result in degrees. Inverse hyperbolic cosine. Inverse cotangent, result in radian. Inverse cotangent, result in degrees. Inverse hyperbolic cotangent. Inverse cosecant, result in radian. Inverse cosecant, result in degrees. 9

ACSCH ASEC ASECD ASECH ASIN ASIND ASINH ATAN ATAN2

Inverse hyperbolic cosecant. Inverse secant, result in radians. Inverse secant, result in degrees. Inverse hyperbolic secant. Inverse sine, result in radians. Inverse sine, result in degrees. Inverse hyperbolic sine. Inverse tangent, result in radians. Four quadrant inverse tangent.

ATAND Inverse tangent, result in degrees. ATANH Inverse hyperbolic tangent. ERFCINV Inverse complementary error function. ERFINV Inverse error function. INV Matrix inverse. PINV Pseudoinverse. IFFT Inverse discrete Fourier transform. IFFT2 Two-dimensional inverse discrete Fourier transform. IFFTN N-dimensional inverse discrete Fourier transform. IFFTSHIFT Inverse FFT shift. inverter.m: %% Inverses of Matrices etc. From this list one can see that the required function is inv. Syntax may then be got from help inv. 1.2.8 The Path Browser

MatLab comes with a large number of M-ﬁles in various directories. 1. When MatLab encounters a name it looks ﬁrst to see if it is a variable name. 2. It then searches for the name as an M-ﬁle in the current directory. (This is one of the reasons to ensure that the program starts in the current directory. 3. It then searches for an M-ﬁle in the directories in the search path. If one of your variables has the same name as an M-ﬁle or a MatLab instruction you will not be able to access that M-ﬁle or MatLab instruction. This is a common cause of problems. The MatLab search path can be added to or changed at any stage by selecting Desktop Tools|Path from the Start Button. Path related functions include addpath Adds a directory to the MatLab search path

10

path Display MatLab search path parh2rc Adds the current path to the MatLab search path rmpath Remove directory from MatLab search path The command cd changes the current working directory 1.2.9 Miscellaneous Commands

Note the following MatLab commands clc Clears the contents of the Command Window clf - Clears the contents of the Figure Window If MATLAB appears to be caught in a loop and is taking too long to ﬁnish a command it may be aborted by ^C (Hold down the Ctrl key and press C). MATLAB will then return to the command prompt diary filename After this command all input and most output is echoed to the speciﬁed ﬁle. The commands diary off and diary on will suspend and resume input to the diary (log) ﬁle.

2

Vectors, Matrices and Arrays

The basic variable in MatLab is an Array. (The numbers entered earlier can be regarded as (1 × 1) arrays, Column vectors as (n × 1) arrays and matrices as (n × m) arrays. MATLAB can also work with multidimensional arrays.

2.1

A Sample MATLAB session

It is recommended that you work through the following sitting at a PC with MATLAB running and enter the commands in the Command window. Most of the calculations involved are simple and they can be checked with a little mental arithmetic. 2.1.1 Entering Matrices

>> x=[1 2 3 4] % assigning values to a (1 by 4) matrix (row vector) x = 1 2 3 4 >> x=[1; 2; 3; 0] % A (4 by 1) (row) vector x = 1

11

2 3 4 >> x=[1,2,3;4,5,6] % (2 by 3) matrix x = 1 4 >> x=[] x = [] 2 5 3 6

%Empty array

%***************************** 2.1.2 Basic Matrix operations

. The following examples are simple. Check the various operations and make sure that you understand them. This will also help you revise some matrix algebra which you will need for your theory. >> x=[1 2;3 4] x = 1 3 >> y=[3 7;5 4] y = 3 5 >> x+y ans = 4 8 >> y-x ans = 9 8 7 4 2 4

%addition of two matrices - same dimensions

%matrix subtraction 2 2 5 0

>> x*y ans =

% matrix multiplication 13 29 15 37

12

Note that when matrices are multiplied their dimensions must conform. The number of columns in the ﬁrst matrix must equal the number of rows in the second otherwise MatLab returns an error. Try the following example. When adding matrices a similar error will be reported if the dimensions do not match >> x=[1 2;3 4] x = 1 2 3 4 >> z=[1,2] z = 1 2 >> x*z ??? Error using ==> mtimes Inner matrix dimensions must agree.

>> inv(x) % find inverse of a matrix ans = -2.00000 1.50000 1.00000 -0.50000

>> x*inv(x) % verify inverse ans = 1.00000 0.00000 0.00000 1.00000

>> y*inv(x) % multiply y by the inverse of x ans = 4.50000 -4.00000 -0.50000 3.00000

>> y/x % alternative expression ans = 4.50000 -4.00000 -0.50000 3.00000

>> inv(x)*y pre-multiply y by the inverse of x ans = 13

1.0e+01 -0.10000 0.20000 >> x\y ans = 1.0e+01 -0.10000 0.20000 2.1.3

* -1.00000 0.85000

% alternative expression - different algorithm - better for regression

* -1.00000 0.85000

Kronecker Product

  a21 B A⊗B = .  .  . an1 B x>> x=[1 2;3 4] x = 1 3 2 4



a11 B

a12 B a22 B . . . an2 B

... ... ...

a1m B a2m B . . . anm B

     

>> I=eye(2,2) I = 1 0 0 1 >> kron(x,I) ans = 1 0 3 0 0 1 0 3 2 0 4 0 0 2 0 4

>> kron(I,x) ans = 14

1 3 0 0 2.1.4

2 4 0 0

0 0 1 3

0 0 2 4

Examples of number formats

>> x=12.345678901234567; >> format loose %includes blank lines to space output >> x x = 12.3457 >> format compact %Suppress blank lines >> x x = 12.3457 >> format long %14 digits after decimal >> x x = 12.34567890123457 >> format short e % exponential or scientific format >> x x = 1.2346e+001 >> format long e >> x x = 1.234567890123457e+001 >> format short g % decimal or exponential >> x x = 12.346 >> format long g >> x x = 12.3456789012346 >> format bank % currency format (2 decimals) >> x

15

x = 12.35 2.1.5 fprintf function

>> fprintf(’%6.2f\n’, x ) 12.35 >> fprintf(’%6.3f\n’, x ) 12.346 >> fprintf(’The number is %6.4f\n’, x ) The number is 12.3457 Here fprintf prints to the command window according to the format speciﬁcation ’%6.4f\n’. In this format speciﬁcation the % indicates the start of a format speciﬁcation. There will be at least 6 digits displayed of which 4 will be decimals in ﬂoating point (f). The \n indicates that the curser will then move to the next line. For more details see page 41. 2.1.6

element by element operations

% .operations >> x=[1 2;3 4]; >> y=[3 7;5 4] >> x .* y %element by element multiplication ans = 3 15 >> y ./ x ans = 3.00000 1.66667 3.50000 1.00000 14 16 %element by element division

>> z=[3 7;0 4]; >> x./z Warning: Divide by zero. ans = 0.3333 0.2857 Inf 1.0000

%mixed scalar matrix operations 16

>> a=2; >> x+a ans = 3 5 >> x-a ans = -1 1 >> x*2 ans = 2 6 >> x/2 ans = 0.50000 1.50000 %exponents % x^a is x^2 or x*x i.e. >> x^a ans = 7 15 10 22h 1.00000 2.00000 4 8 0 2 4 6

% element by element exponent >> z = [1 2;2 1] >> x .^ z ans = 1 9 2.1.7 4 4

miscellaneous functions

Some functions. Operate element by element >> exp(x) 17

ans = 1.0e+01 0.27183 2.00855 >> log(x) ans = 0.00000 1.09861

* 0.73891 5.45982

0.69315 1.38629

>> sqrt(x) ans = 1.00000 1.73205 1.41421 2.00000

Using negative numbers in the argument of logs and square-roots produces an error in many other packages. MATLAB returns complex numbers. Take care!! This is mathematically correct but may not be what you want. >> z=[1 -2] z = 1 -2 >> log(z) ans = 0 >> sqrt(z) ans = 1.0000 >> x-[1 2;3 4] ans = 0 0 0 0 0 + 1.4142i 0.6931 + 3.1416i

>> % determinant >> det(x) ans = -2 >> %trace >> trace(x) ans = 18

5 The function diag(X) where X is a square matrix puts the diagonal of X in a matrix. The function diag(Z) where Z is a column vector outputs a matrix with diagonal Z and zeros elsewhere >> z=diag(x) z = 1 4

>> u=diag(z) u = 1 0 0 4

% Rank of a matrix >> a=[2 4 6 9 3 2 5 4 2 1 7 8] a = 2 3 2 >> rank(a) ans = 3 sum(A) returns sums along diﬀerent dimensions of an array. If A is a vector, sum(A) returns the sum of the elements. If A is a matrix, sum(A) treats the columns of A as vectors, returning a row vector of the sums of each column. >> x=[1 2 3 4] x = 1 >> sum(x) ans = 10 >> sum(x’) ans = 19 2 3 4 4 2 1 6 5 7 9 4 8

10 >> x=[1 2;3 4] x = 1 3 >> sum(x) ans = 4 6 2 4

column-wise from A. An error results if A does not have exactly mn elements >> x=[1 2 3 ; 4 5 6] x = 1 2 3 4 5 6 >> reshape(x,3,2) ans = 1 4 2 5 3 6

The function reshape(A,m,n) returns the m × n matrix B whose elements are taken

blkdiag(A,B,C,) constructs a block diagonal matrix from the matrices A, B c etc. a = 1 3 >> b=5 b = 5 >> c=[6 7 8;9 10 11;12 13 14] c = 6 7 8 9 12 10 13 11 14 2 4

>> blkdiag(a,b,c) ans = 1 3 0 0 0 0 2 4 0 0 0 0 0 0 5 0 0 0 0 0 0 6 9 12 0 0 0 7 10 13 0 0 0 8 11 14 20

This is only a small sample of the available functions eigenvalues and eigenvectors >> A=[54 45 68 A = 54 45 68 45 50 67 68 67 95 45 50 67 68 67 95]

>> eig(A) ans = 0.4109 7.1045 191.4846 >> [V,D]=eig(A) V = 0.3970 0.5898 -0.7032 D = 0.4109 0 0 >> Test=A*V Test = 0.1631 0.2424 -0.2890 >> Test ./ V ans = 0.4109 0.4109 0.4109 >> 7.1045 7.1045 7.1045 191.4846 191.4846 191.4846 0.7631 -0.6378 -0.1042 0 7.1045 0 0.5100 0.4953 0.7033 0 0 191.4846

5.4214 -4.5315 -0.7401

97.6503 94.8336 134.6750

21

2.1.8

sequences

colon operator (:) first:increment:last is a sequence with ﬁrst element ﬁrst second element ﬁrst+ increment and continuing while entry is less than last. >> [1:2:9] ans = 1 >> [2:2:9] ans = 2 >> [1:4] ans = 1 2 3 4 4 6 8

3

5

7

9

>> [1:4]’ ans = 1 2 3 4 %Transpose of a vector >> x x = 1 3 >> x’ ans = 1 2 2.1.9 3 4 2 4

Creating Special Matrices

% creating an Identity Matrix and matrices of ones and zeros >> x=eye(4) 22

x = 1 0 0 0 >> x=ones(4) x = 1 1 1 1 >> x=ones(4,2) x = 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 0 0 0 0 1 0 0 0 0 1

>> x=zeros(3) x = 0 0 0 >> x=zeros(2,3) x = 0 0 >> size(x) ans = 2 2.1.10 3 0 0 0 0 0 0 0 0 0 0

Random number generators

There are two random number generators in standard MATLAB. rand generates uniform random numbers on [0,1) randn generates random numbers from a normal distribution with zero mean and unit 23

variance. >> x=rand(5) x = 0.81551 0.58470 0.70495 0.17658 0.98926 0.55386 0.92263 0.89406 0.44634 0.90159 0.78573 0.78381 0.11670 0.64003 0.52867 0.05959 0.80441 0.45933 0.07634 0.93413 0.61341 0.20930 0.05613 0.14224 0.74421

>> x=rand(5,1) x = 0.21558 0.62703 0.04805 0.20085 0.67641 >> x=randn(1,5) x = 1.29029 2.1.11 1.82176 -0.00236 0.50538 -1.41244

Extracting parts of a matrix, Joining matrices together to get a new larger matrix

>> arr1=[2 4 6 8 10]; >> arr1(3) ans = 6 >> arr2=[1, 2, -3;4, 5, 6;7, 8, 9] arr2 = 1 4 7 >> arr2(2,2) ans = 5 >> arr2(2,:) ans = 4 >> arr2(2,1:2) 24 5 6 2 5 8 -3 6 9

ans = 4 5

>> arr2(end,2:end) ans = 8 2.1.12 9

Using sub-matrices on left hand side of assignment 8 ;9 10 11 12] 3 7 11 4 8 12

>> arr4=[1 2 3 4;5 6 7 arr4 = 1 5 9 2 6 10

>> arr4(1:2,[1,4])=[20,21;22 23] arr4 = 20 22 9 2 6 10 3 7 11 21 23 12

>> arr4=[20,21;22 23] arr4 = 20 22 21 23

>> arr4(1:2,1:2)=1 arr4 = 1 1 1 1 >> arr4=[1 2 3 4;5 6 7 arr4 = 1 5 9 2 6 10 3 7 11 4 8 12

8 ;9 10 11 12]

>> arr4(1:2,1:2)=1 arr4 = 1 1 9 1 1 10 3 7 11 4 8 12

25

2.1.13

Stacking Matrices

>> x=[1 2;3 4] x = 1 3 >> y=[5 6; 7 8] y = 5 7 6 8 2 4

>> z=[x,y,(15:16)’] % join matrices side by side z = 1 3 2 4 5 7 6 8 15 16

>> z=[x’,y’,(15:16)’]’ % Stack matrices vertically z = 1 3 5 7 15 2 4 6 8 16

See also the help ﬁles for the MatLab commands cat, horzcat and vertcat 2.1.14 >> pi pi = 3.1416 >> exp(1) % e = 2.7183 >> clock ans = 1.0e+03 2.00500 YEAR * 0.00100 Month 0.01300 Day 0.02200 hours 0.01100 Minutes 0.01932 Seconds Special Values

2.2

Examples

Work through the following examples using MATLAB. 26

1. let A =

3 5

0 2

and B =

1 4 4 7

Use MATLAB to calculate.

(a) A + B (b) A − B (c) AB (d) AB −1 (e) A/B (f) A\B (g) A. ∗ B (h) A./B (i) A ⊗ B (j) B ⊗ A Use pen, paper and arithmetic to verify that your results are correct.   1 4 3 7    2 6 8 3   Use the MatLab function to show that the rank of  2. Let A =  1 3 4 5    4 13 15 15 A is three. Why is it not four?     14 1 4 3 7        2 6 8 3   and a =  8   3. Solve Ax = a for x where A =   10      1 3 4 5  18 2 1 7 6 trace and determinant of A. Use MATLAB to verify that

4. Generate A which is a 4 × 4 matrix of uniform random numbers. Calculate the (a) The product of the eigenvalues of A is equal to the determinant of A (b) The sum of the eigenvalues of A is equal to the trace of A. (You might ﬁnd the MATLAB functions sum() and prod() helpful - please see the relevant help ﬁles). Do these results hold for an arbitrary matrix A. 5. Let A and B be two 4 × 4 matrices of independent N(0,1) random numbers. If tr(A) is the trace of A. Show that (a) tr(A + B) = tr(A)+tr(B) (b) tr(4A) = 4tr(A) (c) tr(A ) = tr(A) (d) tr(BA)=tr(AB) Which of these results hold for arbitrary matrices? Under what conditions would they hold for non-square matrices? 27

2.3

Regression Example

In this Example I shall use the instructions you have already learned to simulate a set of observations from a linear equation and use the simulated observations to estimate the coeﬃcients in the equation. In the equation yt is related to x2t and x3t according to the following linear relationship.

yt = β1 + β2 x2t + β3 x3t + εt , or in matrix notation y = Xβ + ε where

t = 1, 2, . . . , N

• x2 is a trend variable which takes the values (1,2, . . . 30) • x3 is a random variable with uniform distribution on [3, 5] • εt are independent identically distributed normal random variables with zero mean and constant variance σ 2 . • β1 = 5, β2 = 1 and β3 = 0.1 and εt are iidn(0,.04) (σ 2 = 0.04) 1. Verify that the model may be estimated by OLS. 2. Use MatLab to simulate 50 observations of each of x3 and εt and thus of xt . 3. Using the simulated values ﬁnd OLS estimates of β 4. Estimate the covariance matrix of β and thus the t-statistics for a test that the coeﬃcients of β are zero. 5. Estimate the standard error or the estimate of y 6. Calculate the F-statistic for the signiﬁcance of the regression 7. Export the data to STATA and verify your result. 8. In a simulation exercise such as this two diﬀerent runs will not produce the same result. Any answers submitted should be concise and short and should contain 9. A copy of the m-ﬁle used in the analysis. This should contain comments to explain what is being done 10. A short document giving the results of one simulation and any comments on the results. You might also include the regression table from the STATA analysis. This document should be less than one page in length.

28

A sample answer follows. First the program, then the output and ﬁnally some explanatory notes % example1.m % Regression Example Using Simulated Data %John C Frain %19 November 2006 %values for simulation nsimul=50; beta=[5,1,.1]’; % % Step 1 Prepare and process data for X and y matrices/vectors % x1=ones(nsimul,1); %constant x2=[1:nsimul]’; %trend x3=rand(nsimul,1)*2 +3; X=[x1,x2,x3]; e=randn(nsimul,1)*.2; y= X * beta +e ; % [nobs,nvar]=size(X); % % Estimate Model Note that I have named my estimated variables ols.betahat, ols.yhat, ols.resid etc. The use of the ols. in front of the variable name has two uses. First if I want to do two diﬀerent estimate I will call the estimates ols1. and ols2. or IV. etc. and I can easily put the in a summary table. Secondly this structure has a meaning that is useful in a more advanced use of MATLAB. ols.betahat=(X’*X)\X’*y % Coefficients ols.yhat = X * ols.betahat; % beta(1)*x1+beta(2)*x2+beta(3)*x; ols.resid = y - ols.yhat; % residuals ols.ssr = ols.resid’*ols.resid; % Sum of Squared Residuals ols.sigmasq = ols.ssr/(nobs-nvar); % Estimate of variance ols.covbeta=ols.sigmasq*inv(X’*X); % Covariance of beta ols.sdbeta=sqrt(diag(ols.covbeta));% st. dev of beta ols.tbeta = ols.betahat ./ ols.sdbeta; % t-statistics of beta ym = y - mean(y); ols.rsqr1 = ols.ssr; ols.rsqr2 = ym’*ym; ols.rsqr = 1.0 - ols.rsqr1/ols.rsqr2; % r-squared ols.rsqr1 = ols.rsqr1/(nobs-nvar); 29 % Uniform(3,5)

% N(0,.04) %5*x1 + x2 + .1*x3 + e;

ols.rsqr2 = ols.rsqr2/(nobs-1.0); if ols.rsqr2 ~= 0; ols.rbar = 1 - (ols.rsqr1/ols.rsqr2); % rbar-squared else ols.rbar = ols.rsqr; end; ols.ediff = ols.resid(2:nobs) - ols.resid(1:nobs-1); ols.dw = (ols.ediff’*ols.ediff)/ols.ssr; % durbin-watson fprintf(’R-squared = %9.4f \n’,ols.rsqr); fprintf(’Rbar-squared = %9.4f \n’,ols.rbar); fprintf(’sigma^2 = %9.4f \n’,ols.sigmasq); fprintf(’S.E of estimate= %9.4f \n’,sqrt(ols.sigmasq)); fprintf(’Durbin-Watson fprintf(’Nobs, Nvars = %9.4f \n’,ols.dw); = %6d,%6d \n’,nobs,nvar);

fprintf(’****************************************************\n \n’); % now print coefficient estimates, SE, t-statistics and probabilities %tout = tdis_prb(tbeta,nobs-nvar); % find t-stat probabilities - no %tdis_prb in basic MATLAB - requires leSage toolbox %tmp = [beta sdbeta tbeta tout]; % matrix to be printed tmp = [ols.betahat ols.sdbeta ols.tbeta]; % column labels for printing results namestr = ’ Variable’; bstring = ’ Coef.’; sdstring= ’Std. Err.’; tstring = ’ t-stat.’; cnames = strvcat(namestr,bstring,sdstring, tstring); vname = [’Constant’,’Trend’ ’Variable2’]; The fprintf is used to produce formatted output. See subsection 3.6 fprintf(’%12s %12s %12s %12s \n’,namestr, ... bstring,sdstring,tstring) fprintf(’%12s %12.6f %12.6f %12.6f \n’,... ’ Const’,... ols.betahat(1),ols.sdbeta(1),ols.tbeta(1)) fprintf(’%12s %12.6f %12.6f %12.6f \n’,... ’ Trend’,... ols.betahat(2),ols.sdbeta(2),ols.tbeta(2)) fprintf(’%12s %12.6f %12.6f %12.6f \n’,... ’ Var2’,... ols.betahat(3),ols.sdbeta(3),ols.tbeta(3)) % matrix to be printed

The output of this program should look like 30

R-squared Rbar-squared sigma^2

= = =

0.9998 0.9998 0.0404 0.2010 1.4445

S.E of estimate= Durbin-Watson =

Nobs, Nvars = 50, 3 **************************************************** Variable Const Trend Var2 >> Your answers will of course be diﬀerent Explanatory Notes Most of your MATLAB scripts or programs will consist of three parts 1. Get and Process data Read in your data and prepare vectors or matrices of your left hand side (y), Right hand side (X) and Instrumental Variables (Z) ˆ 2. Estimation Some form of calculation(s) like β = (X X −1 )X y implemented by a MATLAB instruction like betahat = (X’*X)\X*y (where X and y have been set up in the previous step) and estimate of required variances, covariances, standard errors etc. 3. Report Output tables and Graphs in a form suitable for inclusion in a report. 4. Run the program with a smaller number of replications (say 25) and see how the t-statistic on y3 falls. Rerun it with a larger number of replications and see how it rises. Experiment to ﬁnd how many observations are required to get a signiﬁcant coeﬃcient for y3 often. Suggest a use of this kind of analysis. Coef. 4.804620 0.996838 0.147958 Std. Err. 0.229091 0.002070 0.052228 t-stat. 20.972540 481.655756 2.832955

2.4

Simulation – Sample Size and OLS Estimates

This exercise is a study of the eﬀect of sample size on the estimates of the coeﬃcient in an OLS regression. The x values for the regression have been generated as uniform random numbers on the interval [0,100). The residuals are simulated standardised normal random variables. The process is repeated for sample sizes of 20, 100 500 and 2500 simulation is repeated 10,000 times. 31

% example2.m % MATLAB Simulation Example %John C Frain %19 November 2006 % ${ The data files x20.csv, x100.csv, x500.csv and x2500.csv were generated $} using the code below

%Generate Data x20 = 100*rand(20,1) save(’x20.csv’,’x20’,’-ASCII’,’-double’) x100 = 100*rand(100,1) save(’x100.csv’,’x100’,’-ASCII’,’-double’) x500 = 100*rand(500,1) save(’x500.csv’,’x500’,’-ASCII’,’-double’) x2500 = 100*rand(200,1) save(’x2500.csv’,’x2500’,’-ASCII’,’-double’) %} clear nsimul=10000; BETA20=zeros(nsimul,1); % vector - results of simulations with 20 obs. x=load(’-ascii’, ’x20.csv’); % load xdata X=[ones(size(x,1),1),x]; % X matrix note upper case X beta = [ 10;2]; % true values of coefficients % for ii = 1 : nsimul; eps = 20.0*randn(size(X,1),1); % simulated error term y = X * beta + eps; % y values betahat = (X’*X)\X’*y; % estimate of beta BETA20(ii,1)=betahat(2); end fprintf(’Mean and st. dev of 20 obs simulation %6.3f %6.3f\n’ ... ,mean(BETA20),std(BETA20)) %hist(BETA,100) BETA100=zeros(nsimul,1); x=load(’-ascii’, ’x100.csv’); % load xdata X=[ones(size(x,1),1),x]; % X matrix note upper case X beta = [ 10;2]; % % true values of coefficients

for ii = 1 : nsimul;

32

eps = 20.0*randn(size(X,1),1); % simulated error term y = X * beta + eps; % y values betahat = inv(X’*X)*X’*y; % estimate of beta BETA100(ii,1)=betahat(2); end fprintf(’Mean and st. dev of 100 obs simulation %6.3f %6.3f\n’, ... mean(BETA100),std(BETA100)) BETA500=zeros(nsimul,1); x=load(’-ascii’, ’x500.csv’); % load xdata X=[ones(size(x,1),1),x]; % X matrix note upper case X beta = [ 10;2]; % true values of coefficients % for ii = 1 : nsimul; eps = 20.0*randn(size(X,1),1); % simulated error term y = X * beta + eps; % y values betahat = inv(X’*X)*X’*y; % estimate of beta BETA500(ii,1)=betahat(2); end fprintf(’Mean and st. dev of 500 obs simulation %6.3f %6.3f\n’, ... mean(BETA500),std(BETA500)) BETA2500=zeros(nsimul,1); x=load(’-ascii’, ’x2500.csv’); % load xdata note use of lower case x as vector X=[ones(size(x,1),1),x]; % X matrix note upper case X beta = [ 10;2]; % true values of coefficients % for ii = 1 : nsimul; eps = 20.0*randn(size(X,1),1); % simulated error term y = X * beta + eps; % y values betahat = inv(X’*X)*X’*y; % estimate of beta BETA2500(ii,1)=betahat(2); end fprintf(’Mean and st. dev of 2500 obs simulation %6.3f %6.3f\n’, ... mean(BETA2500),std(BETA2500)) n=hist([BETA20,BETA100,BETA500,BETA2500],1.4:0.01:2.6); plot((1.4:0.01:2.6)’,n/nsimul); h = legend(’Sample 20’,’Sample 100’,’Sample 500’,’Sample 2500’);

The output of this program will look like this. On your screen the graph will display coloured lines. 33

Mean and st. dev of 20 obs simulation 2.000 0.165 Mean and st. dev of 100 obs simulation 2.000 0.065 Mean and st. dev of 500 obs simulation 2.000 0.030 Mean and st. dev of 2500 obs simulation
0.14 Sample 20 Sample 100 Sample 500 Sample 2500

1.999

0.049

0.12

0.1

0.08

0.06

0.04

0.02

0

1.4

1.6

1.8

2

2.2

2.4

2.6

2.8

2.5

Example – Macroeconomic Simulation with Matlab

Problem This example is based on the macroeconomic system in Example 10.3 of Shone (2002). There are 10 equations in this economic model. The equations of the system are as follows

ct = 110 + 0.75ydt ydt = yt − taxt taxt = −80 + 0.2yt it = −4rt gt = 330 et = c t + it + g t yt = et−1 mdt = 20 + 0.25yt − 10rt mst = 470 mdt = mst

34

While the aim in Shone (2002) is to examine the system algebraically, here we examine it numerically. Often this may be the only way to solve the system and Matlab is a suitable tool for this work. The model is too simple to be of any particular use in macroeconomics but it does allow one to illustrate the facilities oﬀered by Matlab for this kind of work. Initialise and Describe Variables N = 15 ; % Number of periods for simulation c = NaN * zeros(N,1); %real consumption tax = NaN * zeros(N,1); %real tax yd = NaN * zeros(N,1); %real disposible income i = NaN * zeros(N,1); % real investment g = NaN * zeros(N,1); % real government expenditure e = NaN * zeros(N,1); % real expenditure y = NaN * zeros(N,1); % real income md = NaN * zeros(N,1); % real money demand ms = NaN * zeros(N,1); %real money supply r = NaN * zeros(N,1); % interest rate Simulate g and ms are the policy variables. t=(1:N)’; % time variable g = 330 * ones(N,1); ms = 470 * ones(N,1); y(1) = 2000; The next step is to simulate the model over the required period. In this case this is achieved by a simple reordering of the equations and inverting the money demand equation to give an interest rate equation. In the general case we might need a routine to solve the set of non linear equations or some routine to maximise a utility function. Note that the loop stops one short of the full period and then does the calculations for the ﬁnal period (excluding the income calculation for the period beyond the end of the sample under consideration). for ii = 1:(N-1) tax(ii) = -80 + 0.2 * y(ii); yd(ii) = y(ii) - tax(ii); c(ii) = 110 + 0.75 * yd(ii); md(ii) = ms(ii); r(ii) = (20 + 0.25* y(ii) -md(ii))/10; % inverting money demand 35

i(ii) = 320 -4 * r(ii); e(ii) = c(ii) + i(ii) + g(ii); y(ii+1) = e(ii); end tax(N) = -80 + 0.2 * y(N); yd(N) = y(N) - tax(N); c(N) = 110 + 0.75 * yd(N); md(N) = ms(N); r(N) = (20 + 0.25* y(N) -md(N))/10; i(N) = 320 -4 * r(N); e(N) = c(N) + i(N) + g(N); Now output results and save y for later use. note that the system is in equilibrium. Note that in printing we use the transpose of base base = [t,y,yd,c,g-tax,i,r]; fprintf(’ t y yd

c

g-tax

i

r\n’)

fprintf(’%7.0f%7.0f%7.0f%7.0f%7.0f%7.0f%7.0f\n’,base’) ybase = y; t 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 y 2000 2000 2000 2000 2000 2000 2000 2000 2000 2000 2000 2000 2000 2000 2000 yd 1680 1680 1680 1680 1680 1680 1680 1680 1680 1680 1680 1680 1680 1680 1680 c 1370 1370 1370 1370 1370 1370 1370 1370 1370 1370 1370 1370 1370 1370 1370 g-tax 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 i 300 300 300 300 300 300 300 300 300 300 300 300 300 300 300 r 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5

Revised Simulation We increase g to 350 and examine the passage to the new equilibrium. Basically we run the same program with a diﬀerent starting value for g.

36

N = 15 ; % Number of periods for simulation c = NaN * zeros(N,1); %real consumption tax = NaN * zeros(N,1); %real tax yd = NaN * zeros(N,1); %real disposible income i = NaN * zeros(N,1); % real investment g = NaN * zeros(N,1); % real government expenditure e = NaN * zeros(N,1); % real expenditure y = NaN * zeros(N,1); % real income md = NaN * zeros(N,1); % real money demand ms = NaN * zeros(N,1); %real money supply r = NaN * zeros(N,1); % interest rate % Policy Variables g = 350 * ones(N,1); ms = 470 * ones(N,1); t=(1:N)’; y(1) = 2000; for ii = 1:(N-1) tax(ii) = -80 + 0.2 * y(ii); yd(ii) = y(ii) - tax(ii); c(ii) = 110 + 0.75 * yd(ii); md(ii) = ms(ii); r(ii) = (20 + 0.25* y(ii) -md(ii))/10; % inverting money demand i(ii) = 320 -4 * r(ii); e(ii) = c(ii) + i(ii) + g(ii); y(ii+1) = e(ii); end tax(N) = -80 + 0.2 * y(N); yd(N) = y(N) - tax(N); c(N) = 110 + 0.75 * yd(N); md(N) = ms(N); r(N) = (20 + 0.25* y(N) -md(N))/10; i(N) = 320 -4 * r(N); e(N) = c(N) + i(N) + g(N);

policy = [t,y,yd,c,g-tax,i,r]; fprintf(’ t y yd c g-tax i r\n’) fprintf(’%7.0f%7.0f%7.0f%7.0f%7.0f%7.0f%7.0f\n’,policy’) ypolicy =y;

37

t 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

y 2000 2020 2030 2035 2038 2039 2039 2040 2040 2040 2040 2040 2040 2040 2040

yd 1680 1696 1704 1708 1710 1711 1712 1712 1712 1712 1712 1712 1712 1712 1712

c 1370 1382 1388 1391 1393 1393 1394 1394 1394 1394 1394 1394 1394 1394 1394

g-tax 30 26 24 23 23 22 22 22 22 22 22 22 22 22 22

i 300 298 297 297 296 296 296 296 296 296 296 296 296 296 296

r 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6

Now we compare results in a table and a graph. Note that income converges to a new limit. fprintf(’ t ybase ypolicy\n’) fprintf(’%7.0f%7.0f%7.0f\n’,[t,ybase, ypolicy]’) plot(t,[ybase,ypolicy]) title(’Equilibrium and Shocked Macro-system’) xlabel(’Period’) ylabel(’Income’) axis([1 15 1995 2045]) t 1 2 3 4 5 6 7 8 9 10 11 12 13 ybase ypolicy 2000 2000 2000 2000 2000 2000 2000 2000 2000 2000 2000 2000 2000 2000 2020 2030 2035 2038 2039 2039 2040 2040 2040 2040 2040 2040 38

14 15

2000 2000

2040 2040
Equilibrium and Shocked Macro−system

2045 2040 2035 2030 2025 Income 2020 2015 2010 2005 2000 1995 2 4 6 8 Period 10 12 14

3
3.1

Data input/output
Native MatLab data ﬁles

The instruction Save filename saves the contents of the workspace in the ﬁle ’filename.mat’. save used in the default manner saves the data in a binary format. The instruction save filename, var1, var2 saves var1 and var2 in the ﬁle filename.mat. Similarly the commands Load filename and load filename, var1, var2. load the contents of ’filename.mat’ or the speciﬁed variables from the ﬁle into the workspace. In general .mat ﬁles are nor easily readable in most other packages. They are ideal for use within MATLAB and for exchange between MATLAB users. (note that there may be some incompatibilities between diﬀerent versions of MATLAB). These .mat ﬁles are binary and can not be examined in a text editor. .mat is the default extension for a MATLAB data ﬁle. If you use another extension, say .ext the option Save mat filename.ext should be used with the save and load commands. It is possible to use save and load to save and load text ﬁles but these instructions are very limited. If your data are in EXCEL or csv format the methods described below are better

39

3.2

Importing from Excel

The sample ﬁle g10xrate.xls contains daily observations on the exchange rates of G10 countries and we wish to analyse them with MATLAB. There are 6237 observations of each exchange rate in the columns of the EXCEL ﬁle. The easiest way to import these data into MATLAB is to use the File|import data wizard and follow the prompts. In this case the import wizard did not pick out the series names from the Excel ﬁle. I imported the entire data matrix as a matrix and extracted the individual series in MATLAB. One can save the data in MATLAB format for future use in MATLAB This is illustrated in the code below. USXJPN = data(:,1); USXFRA = data(:,2); USXSUI = data(:,3); USXNLD = data(:,4); USXUK = data(:,5); USXBEL = data(:,6); USXGER = data(:,7); USXSWE = data(:,8); USXCAN = data(:,9); USXITA = data(:,10); save(’g10xrate’, ’USXJPN’,’USXFRA’,’USXSUI’,’USXNLD’,’USXUK’,’USXBEL’,’USXGER’, ... ’USXSWE’,’USXCAN’,’USXITA’) Note that I have listed the series to be saved as I did not wish to save the data matrix. The same eﬀect could have been achieved with the uiimport command.

3.3

Reading from text ﬁles

The import wizard can also import many types of text ﬁle including the comma separated ﬁles we have used in STATA. The missing data code in Excel csv ﬁles is #NA. The version of MATLAB that i am using has problems reading this missing value code and it should be changed to NaN (the MATLAB missing value code) before importing csv data. In this case the import wizard recognised the column names. It is important that you check that all your data has been imported correctly. The MATLAB functions textscan or textread can read various text ﬁles and allow a greater degree of ﬂexibility than that available from uiimport. This ﬂexibility is obtained at a cost of greater complexity. Details are given in the Help ﬁles. I would think that most users will not need this ﬂexibility but it is there if needed.

40

3.4

Exporting data to EXCEL, STATA and other programs

The command xlswrite(’filename’,M) writes the matrix M to the ﬁle filename in the current working directory. If M is n × m the numeric values in the matrix are written to the ﬁrst n row and m columns in the ﬁrst sheet in the spreadsheet. The command csvwrite(’filename’,M) writes the matrix M to the ﬁle filename in the current working directory. You can use this ﬁle to transfer data to STATA. Alternatively export your Excel ﬁle from Excel in csv format.

3.5

Stat/Transfer

Another alternative is to use the Stat/Transfer package which allows the transfer of data ﬁles between a large number of statistical packages.

3.6

Formatted Output

The MATLAB function fprintf() may be used to produce formatted output on screen1 . The following MATLAB program gives an example of the us of the fprintf() function. Sample MATLAB program demonstrating Formatted Output clear degrees_c =10:10:100; degrees_f = (degrees_c * 9 /5) + 32; fprintf(’\n\n Conversion from degrees Celsius \n’); fprintf(’ fprintf(’ to degrees Fahrenheit\n\n’ ); Celsius Fahrenheit\n’);

for ii = 1:10; fprintf(’%12.2f%12.2f\n’,degrees_c(ii),degrees_f(ii)); end % fprintf(... ’\n\n%5.2f degrees Celsius is equivalent of %5.3f degrees fahrenheit\n’, ... degrees_c(1),degrees_f(1)) Output of Sample MATLAB program demonstrating Formatted Output

Conversion from degrees Celsius is only one of a large number of C-style input/output functions in C. These allow considerable ﬂexibility in sending formatted material to the screen of to a ﬁle. The MATLAB help ﬁles give details of the facilities available. If further information is required one might consult a standard test book on the C programming language
1 fprintf()

41

to degrees Fahrenheit Celsius 10.00 20.00 30.00 40.00 50.00 60.00 70.00 80.00 90.00 100.00 Fahrenheit 50.00 68.00 86.00 104.00 122.00 140.00 158.00 176.00 194.00 212.00

10.00 degrees Celsius is equivalent of 50.000 degrees fahrenheit Note the following • The ﬁrst argument of the fprintf() function is a kind of format statement included within ’ marks. • The remaining arguments are a list of variables or items to be printed separated by commas • Within the format string there is text which is produced exactly as set down. There are also statements like %m.nf which produces a decimal number which is allowed m columns of output and has n places of decimals. These are applied in turn to the items in the list to be printed. • This f format is used to output ﬂoating point numbers there are a considerable number or other speciﬁers to output characters, strings, and number in formats other than ﬂoating point. • If the list to be printed is too long the formats are recycled. • Not the use of \n which means skip to the next line. This is essential.

3.7

Producing material for inclusion in a paper

A considerable amount of the material in this note was produced from MATLAB mﬁles using the —File—Publish to— facilities in the MATLAB m-ﬁle editor which A produces output in WORD, Powerpoint, L TEX , HTML etc. for inclusion in papers, presentations etc. The facilities are described in the help ﬁles and may vary from version to version of Matlab. 42

To Access these facilities you must ﬁrst turn them on in the Matlab Editor. This is done by —Cell—Enable Cell Mode— in the editor menu. Cell mode enables you to divide your m-ﬁle into cells or sections. (Do not confuse cell mode in the editor with cell data structures in Matlab. It is unfortunate that these two diﬀerent concepts have the same name) Enabling cell mode adds a new set of buttons to the menu bar and enables a set of items in the cell menu item. The menu items allow one to • Disable cell mode • Evaluate current cell • Evaluate current cell and advance • Evaluate entire ﬁle • Insert Cell Divider • Insert Cell Dividers around Selection • Insert Cell Markup – Cell Title – Descriptive Text – Bold Text – Monospaced Text – Preformatted Text – Bulleted Text – TEX Equation • Next Cell • Previous Cell Section 2.5, including its graph, was subject to very minor editing before being added to this document. I have also used the facility to produce transparencies for lectures on MATLAB.

4

Decision and Loop Structures.

There are four basic control (Decision or Loop Structures) available in MATLAB if statements The basic form of the if statement is if conditions statements end

43

The statements are only processed if the conditions are true The conditions can include the following operators == ∼= equal not equal less than greater than less than or equal to greater than or equal to logical and logical and (for scalars) short-circuiting logical or logical or and (for scalars) short-circuiting logical exclusive or true if all elements of vector are nonzero true if any element of vector is nonzero

< >

= & && | || xor

all any

The if statement may be extended if conditions statements1 else statements2 end

in which case statements1 are used if conditions are true and statements2 if false. This if statement may be extended again if conditions1 statements1 elseif conditions2 statements2 else statements3 end

with an obvious meaning (I hope). for The basic form of the for group is for variable = expression statements end 44

Here expression is probably a vector. statements is processed for each of the values in expression. The following example shows the use of a loop within a loop >> for ii = 1:3 for jj=1:3 total=ii+jj; fprintf(’%d + %d = %d \n’,ii,jj,total) end end 1 + 1 = 2 1 + 2 = 3 1 + 3 = 4 2 + 1 = 3 2 + 2 = 4 2 + 3 = 5 3 + 1 = 4 3 + 2 = 5 3 + 3 = 6 while The format of the while statement is while conditions statements end The while statement has the same basic functionality as the for statement. The for statement will be used when one knows precisely when and how many times an operation will be repeated. The statements are repeated so long as conditions are true switch An example of the use of the switch statement follows switch p case 1 x = 24 case 2 x = 19 case 3 x = 15 otherwise error(’p must be 1, 2 or 3’) end Use matrix statements in preference to loops. Not only are they more eﬃcient but they are generally easier to use. That said there are occasions where one can not use a matrix statement. 45

If you wish to ﬁll the elements of a vector or matrix using a loop it is good practice to initialise the vector or matrix ﬁrst. For example if you wish to ﬁll a 100 × 20 matrix, X, using a series of loops one could initialise the matrix using one of the following commands X = ones(100,20) X = zeros(100,20) X = ones(100,20)*NaN X = NaN(100,20)

5

Elementary Plots

Simple graphs can be produced easily in MatLab. The following sequence %values for simulation nsimul=50; beta=[5,1,.1]’; % x1=ones(nsimul,1); %constant x2=[1:nsimul]’; %trend x3=rand(nsimul,1)*2 +3; % Uniform(3,5) x=[x1,x2,x3]; e=randn(nsimul,1)*.2; y= x * beta +e ; % [nobs,nvar]=size(x); betahat=inv(x’*x)*x’*y %g yhat = x * betahat; % beta(1)*x1-beta(2)*x2-beta(3)*x; resid = y - yhat; plot(x2,resid) title(’Graph Title’) xlabel(’Time’) ylabel(’Residual’) % N(0,.04) %5*x1 + x2 + .1*x3 + e;

46

Graph Title 0.5 0.4 0.3 0.2 0.1 Residual 0 −0.1 −0.2 −0.3 −0.4 −0.5 0 10 20 Time repeats the earlier OLS simulation, opens a graph window, draws a graph of the residuals against the trend in the ols-simulation exercise, puts a title on the graph and labels the x and y axes. The vectors x2 and resid must have the same dimensions. This graph was saved in eps format and imported into this document.

30

40

50

6
6.1

Systems of Regression Equations
Using Matlab to estimate systems of regression equations

This section contains two examples of the estimation of systems of equations. The ﬁrst is an examination of the classic Grunfeld investment data set. Many textbooks use this dataset to illustrate various features of system estimation. Green (2000) is the source of the data used here. Later editions of this book also examine these data but in less detail. The MATLAB output also includes corresponding analysis using the le Sage econometrics package which is covered in section 8 of this note. As an exercise the user might extend the analysis to include various Likelihood Ratio tests of the restrictions imposed by the various estimation procedures.

47

Analysis of Grunfeld Investment data Introduction The basic system model that we are looking at here is yti = Xti βi + εti where 1 ≤ i ≤ M represents the individual agent or country or item for which we are estimating some equation and 1 ≤ t ≤ T represents the tth measurement on the ith unit.
2 We assume that the variance of εti , σi is constant for 1 ≤ t ≤ T . Each Xi is T × ki . We may write these equations

y1 = X1 β1 + ε1 y2 = X2 β2 + ε2 ... yM = XM βM + εM In this section we will assume that X is exogenous. By imposing various cross-equation restrictions on the βi and the covariances of the εti we obtain a variety of estimators (e.g. Pooled OLS, Equation be Equation OLS, SUR). The variables included in the Grunfeld analysis are • • • • • FIRM : There are 10 ﬁrms YEAR : Data are from 1935 to 1954 (20 years) I : Gross Investment F : Value of Firm C : Stock of Plant and Equipment

For more details see Green (2000, 2008) or the original references listed there. • • • • • • To reduce the amount of detail we shall restrict analysis to 5 ﬁrms Firm no 1 : GM - General Motors Firm no 4 : GE - general electric Firm no 3 : CH - Chrysler Firm no 8 : WE - Westinghouse Firm no 2 : US - US Steel

test names on the data set are not imported as I wish to deﬁne these myself. This sets up a matrix containing data. I save data in Matlab form the ﬁrst time. Later I use load data; % to reload the data as below 48

To start the analysis is use the MATLAB Import data using |File|Import Data]. The

Load and Generate data load data Y_GM = data(1:20, 3); % I X_GM = [ones(20,1),data(1:20,[4 5])]; % constant F C Y_GE = data(61:80, 3); % I X_GE = [ones(20,1),data(61:80,[4 5])]; % constant F C Y_CH = data(41:60, 3); % I X_CH = [ones(20,1),data(41:60,[4 5])]; % constant F C Y_WE = data(141:160, 3); % I X_WE = [ones(20,1),data(141:160,[4 5])]; % constant F C Y_US = data(21:40, 3); % I X_US = [ones(20,1),data(21:40,[4 5])]; % constant F C We now estimate the coeﬃcients imposing various restrictions. Each estimation involves the following steps 1. Set up the required y and X matrices. 2. Estimate the required coeﬃcients. 3. Estimate standard errors, t-statistics etc. 4. Report. Pooled OLS The restrictions imposed by Pooled OLS are that corresponding coeﬃcients are the same across equations. We also assume that the variance of the disturbances is constant across equations.2 Thus, in this case ki = k, for all i We can therefore assume that each observation on each unit is one more observation from the same single equation system. We may write the system as follows     X1 y1      y2   X2   =  β  ...   ...      XM yM (M T × 1) (M T × k)  ε1    ε2     ...    εM 

+

(k × 1)

(M T × 1)

or, more compactly, using the obvious notation y = Xβ + ε ˆ and β may be estimated by β = (X X)−1 X y etc. This is implemented in MATLAB as follows –
2 We

could of course relax this condition and estimate Heteroskedastic Consistent Standard Errors

49

Y = [Y_GM’, Y_GE’, Y_CH’, Y_WE’, Y_US’]’; % leave out ; for testing if % necessary delete during run or output will be unreadable X = [X_GM’, X_GE’, X_CH’, X_WE’, X_US’]’; pols.beta = (X’*X)\X’*Y; pols.uhat = Y - X*pols.beta ; pols.sigsq = (pols.uhat’*pols.uhat)/(size(X,1)-size(X,2));%(T-k) pols.sdbeta = sqrt(diag(inv(X’*X))*pols.sigsq); pols.tbeta = pols.beta ./ pols.sdbeta; pols.se = sqrt(pols.sigsq); label = [’Constant ’; ’F ’; ’C disp(’OLS Results using stacked matrices’) disp(’ for ii=1:size(X,2) coef sd t-stat’) ’];

fprintf(’%s%10.4f%10.4f%10.4f\n’,label(ii,:),pols.beta(ii),pols.sdbeta(ii), pols.tbeta(ii)) end fprintf(’Estimated Standard Error %10.4f\n\n\n’,pols.se) OLS Results using stacked matrices Constant F C coef -47.0237 0.1059 0.3014 sd 21.6292 0.0114 0.0437 t-stat -2.1741 9.2497 6.8915

Estimated Standard Error %

128.1429

% Verification using Lesage package % pooled = ols(Y, X); vnames= [’I ’Constant ’F ’C prt(pooled,vnames) ’; ’; ’; ’];

\begin{verbatim} Ordinary Least-squares Estimates Dependent Variable = I R-squared = 0.7762 Rbar-squared sigma^2 Durbin-Watson = 0.7716 = 16420.6075 = 0.3533

50

Nobs, Nvars = 100, 3 *************************************************************** Variable Coefficient t-statistic t-probability Constant F C -47.023691 0.105885 0.301385 -2.174080 9.249659 6.891475 0.032132 0.000000 0.000000

Equation by equation OLS This section assumes that the coeﬃcients vary across units. (In panel data estimation we assume that only the constant terms vary across units). We also assume that there is no contemporaneous correlation between the disturbances in the system. We may write the system of equations as

or more compactly using the obvious substitutions y = Xβ + ε

  X1 y1     y2   0  . = .  .   .  .   . 0 yM 

0 X2 . . . 0

··· ··· .. . ···

   ε1 0    0   ε2  . β +  .   .   .   .  . εM XM

where y and ε are T M × 1, X is T M × kM and β is kM × 1. y, ε, and β are stacked versions of yi , εi , and βi . The variance of ε is given by Ω = E [εε ]  2 σ 1 IT   0 = .  .  . 0
2 σ1   0 =   0

0
2 σ 2 IT . . .

··· .. .

···

0 0 . . .
2 σ M IT 

     

0 0
2 σ2

··· 0



0 0

··· .. .

···

0

···

2 σM

 0   ⊗ IT  0 

The coding of this example should be relatively clear. Perhaps the most diﬃcult part is the estimation of the variances. The procedure here is very similar to that used in the ﬁrst step of the SUR estimation procedure except that the contemporaneous correlation is used to improve the estimates. It should be noted that, in this case diﬀerent variables are likely to be used in diﬀerent equations, The only change required is in the calculation of the X matrix. 51

% Y as before X=blkdiag(X_GM ,X_GE , X_CH , X_WE , X_US); eqols.beta = (X’*X)\X’*Y; eqols.uhat = Y - X*eqols.beta ; eqols.temp = reshape(eqols.uhat,size(X_GM,1),5); %residuals for % each firm in a column eqols.sigsq1 =eqols.temp’*eqols.temp/(size(X_GM,1)-size(X_GM,2)); eqols.sigsq = diag(diag(eqols.sigsq1)); % Remove non-diagonal elements %eqols.sdbeta = sqrt(diag(inv(X’*X)*X’*kron(eye(size(X_GM,1)),eqols.sigsq)*X*inv(X’*X))); eqols.covarbeta = inv(X’*X)*kron(eqols.sigsq,eye(3)); eqols.sdbeta = diag(sqrt(eqols.covarbeta)); eqols.tbeta=eqols.beta ./ eqols.sdbeta; eqols.se=sqrt(diag(eqols.sigsq)); % % Write results % disp(’OLS equation by equation using stacked matrices’) disp(’OLS estimates GE equation’) firm = [’GE’; ’GM’; ’CH’; ’wE’; ’US’]; for jj = 1:5 % Loop over firms fprintf(’\n\n\n’) disp(’ coef sd for ii=1:3 %Loop over coefficients t-stat’)

fprintf(’%10s%10.4f%10.4f%10.4f\n’,label(ii), ... eqols.beta(ii+(jj-1)*3),eqols.sdbeta(ii+(jj-1)*3), ... eqols.tbeta(ii+(jj-1)*3)) end fprintf(’Standard Error is %10.4f\n’,eqols.se(jj)); end OLS equation by equation using stacked matrices OLS estimates GE equation

coef C -149.7825 F 0.1193

sd 105.8421 0.0258

t-stat -1.4151 4.6172 52

C 0.3714 0.0371 Standard Error is 91.7817

10.0193

C F C

coef -6.1900 0.0779 0.3157

sd 13.5065 0.0200 0.0288 13.2786

t-stat -0.4583 3.9026 10.9574

Standard Error is

coef C F -9.9563 0.0266

sd 31.3742 0.0156

t-stat -0.3173 1.7057 5.9015

C 0.1517 0.0257 Standard Error is 27.8827

C F C

coef -0.5094 0.0529 0.0924

sd 8.0153 0.0157 0.0561 10.2131

t-stat -0.0636 3.3677 1.6472

Standard Error is

coef C F -49.1983 0.1749

sd 148.0754 0.0742

t-stat -0.3323 2.3566 2.7369

C 0.3896 0.1424 Standard Error is 96.4345 Verify using le Sage Toolbox olsestim=ols(Y_US,X_US); prt(olsestim, vnames); Ordinary Least-squares Estimates Dependent Variable = I R-squared = 0.4709

53

Rbar-squared sigma^2 Durbin-Watson

= 0.4086 = 9299.6040 = 0.9456

Nobs, Nvars = 20, 3 *************************************************************** Variable Constant F C Coefficient -49.198322 0.174856 0.389642 t-statistic -0.332252 2.356612 2.736886 t-probability 0.743761 0.030699 0.014049

SUR Estimates
Suppose that we have a random sample of households and we are have time series data on expenditure on holidays (yit ) and relevant explanatory variables. Suppose that we have suﬃcient data to estimate a single equation for each person in the sample. We also assume that there is no autocorrelation in each equation (often a rather heroic assumption). During the peak of the business cycle it is likely that many of the persons in the sample spend above what they do at the trough. Thus it is likely that there will be contemporaneous correlation between the errors in the system.  σ if i = j ij E[εti εsj ] = 0 if i = j  σ11 σ12 σ22 . . . σM 2 ··· σ1M 

Thus we may write the contemporaneous covariance matrix (Σ) as

  σ21 Σ= .  .  .  Σ  0 Ω=. . . 0

··· .. .

σM 1

···

 σ2M  .  .  .  σM M

The total covariance matrix is, in this case, given by 0 Σ . . . 0 ··· ··· .. . ···  0  0 .  .  .  Σ

= Σ ⊗ IT

If Ω were known we would use GLS to get optimum estimates of β. In this case we can obtain a consistent estimate of Σ from the residuals in the equation by equation OLS estimate that we have just completed. We can then use this consistent estimate in Feasible GLS.

54

Omega = kron(eqols.sigsq1,eye(20,20)); % Page page 256 eqsur.beta= inv(X’*inv(Omega)*X)*X’*inv(Omega)*Y; eqsur.yhat = X * eqsur.beta; eqsur.uhat = Y - eqsur.yhat; eqsur.temp = reshape(eqsur.uhat,20,5); eqsur.omega = eqsur.temp’ * eqsur.temp /size(X_GM,1); %(size(X_GM,1)-size(X_GM,2)); eqsur.covar = inv(X’*inv(kron(eqsur.omega, eye(20)))*X); eqsur.sdbeta = sqrt(diag(eqsur.covar)); eqsur.tbeta = eqsur.beta ./ eqsur.sdbeta; eqsur.se = sqrt(diag(eqsur.omega)); %print results fprintf(’SUR estimates\n’); for jj = 1:5 % Loop over firms fprintf(’\n\n\n’) disp(’ coef sd for ii=1:3 %Loop over coefficients t-stat’)

fprintf(’%s%10.4f%10.4f%10.4f\n’,label(ii), ... eqsur.beta(ii+(jj-1)*3),eqsur.sdbeta(ii+(jj-1)*3), ... eqsur.tbeta(ii+(jj-1)*3)) end fprintf(’Standard Error is %10.4f\n’,eqsur.se(jj)); end SUR estimates

coef C -168.1134 F 0.1219 84.9017 0.0204

sd -1.9801 5.9700

t-stat

C 0.3822 0.0321 11.9109 Standard Error is 84.9836

C

0.9980

coef 11.5661

sd 0.0863

t-stat

F 0.0689 0.0170 4.0473 C 0.3084 0.0260 11.8766 Standard Error is 12.3789

55

C F

-21.1374 0.0371

coef 24.3479 0.0115

sd -0.8681 3.2327

t-stat

C 0.1287 0.0212 6.0728 Standard Error is 26.5467

coef C F C 1.4075 0.0564 0.0429 5.9944 0.0106 0.0382

sd 0.2348 5.3193 1.1233 9.7420

t-stat

Standard Error is

C F

62.2563 0.1214

coef 93.0441 0.0451

sd 0.6691 2.6948

t-stat

C 0.3691 0.1070 3.4494 Standard Error is 90.4117

SUR in LeSage toolbox y(1).eq = Y_GM; y(2).eq = Y_GE; y(3).eq = Y_CH; y(4).eq = Y_WE; y(5).eq = Y_US; XX(1).eq = X_GM; XX(2).eq = X_GE; XX(3).eq = X_CH; XX(4).eq = X_WE; XX(5).eq = X_US; neqs=5; sur_result=sur(neqs,y,XX); prt(sur_result) Seemingly Unrelated Regression -- Equation System R-sqr R-squared Rbar-squared = = = 0.8694 0.9207 0.9113 1

56

sigma^2 Durbin-Watson Nobs, Nvars

= 458183.2995 = 0.0400 = 20, 3

*************************************************************** Variable Coefficient t-statistic t-probability variable variable variable 1 2 3 -168.113426 0.121906 0.382167 -1.980094 5.969973 11.910936 0.064116 0.000015 0.000000

Seemingly Unrelated Regression -- Equation System R-sqr = 0.8694 R-squared Rbar-squared sigma^2 Durbin-Watson = = 0.9116 0.9012

2

= 8879.1368 = 0.0310

Nobs, Nvars = 20, 3 *************************************************************** Variable Coefficient t-statistic t-probability variable variable variable 1 2 3 0.997999 0.068861 0.308388 0.086286 4.047270 11.876603 0.932247 0.000837 0.000000

Seemingly Unrelated Regression -- Equation System R-sqr R-squared Rbar-squared sigma^2 Durbin-Watson = = = 0.8694 0.6857 0.6488

3

= 11785.8684 = 0.0202

Nobs, Nvars = 20, 3 *************************************************************** Variable variable variable variable 1 2 3 Coefficient -21.137397 0.037053 0.128687 t-statistic -0.868140 3.232726 6.072805 t-probability 0.397408 0.004891 0.000012

Seemingly Unrelated Regression -- Equation System R-sqr R-squared Rbar-squared = = = 0.8694 0.7264 0.6943

4

57

sigma^2 Durbin-Watson Nobs, Nvars

= 2042.8631 = 0.0323 = 20,

3

*************************************************************** Variable Coefficient t-statistic t-probability variable variable variable 1 2 3 1.407487 0.056356 0.042902 0.234802 5.319333 1.123296 0.817168 0.000056 0.276925

Seemingly Unrelated Regression -- Equation System R-sqr = 0.8694 R-squared Rbar-squared sigma^2 Durbin-Watson = = 0.4528 0.3884

5

= 173504.8346 = 0.0103

Nobs, Nvars = 20, 3 *************************************************************** Variable Coefficient t-statistic t-probability variable variable variable 1 2 3 62.256312 0.121402 0.369111 0.669105 2.694815 3.449403 0.512413 0.015340 0.003062

Cross-equation sig(i,j) estimates equation eq 1 eq 2 eq 3 eq 1 eq 2 eq 3 eq 4 eq 5 7222.2204 -315.6107 -315.6107 153.2369 601.6316 3.1478 129.7644 -2446.3171 601.6316 3.1478 704.7290

eq 4

eq 5

129.7644 -2446.3171 16.6475 414.5298 201.4385 1298.6953 94.9067 613.9925 613.9925 8174.2798

16.6475 201.4385 414.5298 1298.6953

Cross-equation correlations equation eq 1 eq 2 eq 1 eq 2 eq 3 eq 4 eq 5 1.0000 -0.3000 0.2667 0.1567 -0.3184 -0.3000 1.0000 0.0096 0.1380 0.3704

eq 3 0.2667 0.0096 1.0000 0.7789 0.5411

eq 4 0.1567 0.1380 0.7789 1.0000 0.6971

eq 5 -0.3184 0.3704 0.5411 0.6971 1.0000

58

6.2

Exercise – Using Matlab to estimate a simultaneous equation systems

Consider the demand-supply model qt qt = = β11 + β21 xt2 + β31 xt2 + γ21 pt + ut1 β12 + β42 xt4 + β52 xt5 + γ22 pt + ut2 , (1) (2)

where qt is the log of quantity, pt is the log of price, xt2 is the log of income, xt3 is a dummy variable that accounts for demand shifts xt4 and xt5 are input prices. Thus equations (1) and (2) are demand and supply functions respectively. 120 observations generated by this model are in the ﬁle demand-supply.csv 1. Comment on the identiﬁcation of the system. Why can the system not be estimated using equation by equation OLS. For each of the estimates below produce estimates, standard errors and t-statistics of each coeﬃcient. Also produce standard errors for each equation. 2. Estimate the system equation by equation using OLS. 3. Estimate the system equation by equation using 2SLS. Compare the results with the OLS estimates. 4. Set up the matrices of included variables, exogenous variables required to do system estimation. 5. Do OLS estimation using the stacked system and compare results with the equation by equation estimates. 6. Do 2SLS estimation using the stacked system and compare results with the equation by equation estimates. 7. Do 3SLS estimation using the stacked system and compare results with the 2SLS estimates. 8. Comment on the identiﬁcation of the system. 9. How can the method be generalised to estimate other GMM estimators? Estimate the optimum GMM estimator for the system and compare your results with the previous estimators.

7

User written functions in MATLAB

One of the most useful facilities in MATLAB is the facility to write ones own functions and use them in the same way as a native MATLAB functions. We are already familiar with m-ﬁles which contain lists of MATLAB instructions. Such ﬁles are known as script 59

ﬁles and allow us to do repeat an analysis without having to retype all the instructions. Suppose we wish to write a function that estimates the density function of a normal distribution, √ (x − µ)2 1 exp − 2σ 2 2π σ

, we have the following on a ﬁle normdensity.m function f = normdensity(z, mu, sigma); % Calculates the Density Function of the Normal Distribution % with mean mu % and standard deviation sigma % at a point z % sigma must be a positive non-zero real number if sigma > help ols PURPOSE: least-squares regression --------------------------------------------------USAGE: results = ols(y,x) where: y = dependent variable vector (nobs x 1) x = independent variables matrix (nobs x nvar) --------------------------------------------------RETURNS: a structure results.meth = ’ols’ results.beta = bhat results.tstat = t-stats results.bstd results.yhat (nvar x 1) (nvar x 1)

= std deviations for bhat (nvar x 1) = yhat (nobs x 1)

results.resid = residuals (nobs x 1) results.sige = e’*e/(n-k) scalar results.rsqr results.rbar results.dw = rsquared scalar = rbar-squared scalar = Durbin-Watson Statistic 62

results.nobs results.nvar results.y

= nobs = nvars = y data vector (nobs x 1)

results.bint = (nvar x 2 ) vector with 95% confidence intervals on beta --------------------------------------------------SEE ALSO: prt(results), plt(results) --------------------------------------------------Overloaded functions or methods (ones with the same name in other directories) help localmod/ols.m After running the ols function a structure containing the results is available. The variable results.beta contains the estimated β-coeﬃcients, results.tstat their tstatistics, results.bint the 95% conﬁdence intervals for the estimates and similar for the other variables deﬁned. Each estimation command produces its results in a similar structure. To see how to print a summary of these results >> help prt_reg PURPOSE: Prints output using regression results structures --------------------------------------------------USAGE: prt_reg(results,vnames,fid) Where: results = a structure returned by a regression vnames = an optional vector of variable names = optional file-id for printing results to a file (defaults to the MATLAB command window) --------------------------------------------------NOTES: e.g. vnames = strvcat(’y’,’const’,’x1’,’x2’); e.g. fid = fopen(’ols.out’,’wr’); use prt_reg(results,[],fid) to print to a file with no vnames -------------------------------------------------RETURNS: nothing, just prints the regression results -------------------------------------------------SEE ALSO: prt, plt --------------------------------------------------Thus to display the results of the previous regression on the screen in MATLAB one would enter3 prt_reg(result) in this context is the name of a MATLAB variable and one could substitute for result any valid MATLAB variable name.
3 result

fid

63

Availability on other PCs The LeSage toolbox is available for download free on the internet. The only requirement on the package web site is that “Anyone is free to use these routines, no attribution (or blame) need be placed on the author/authors.” The econometrics package is not available by default when MATLAB is installed on a PC. It may be downloaded from http://www.spatial-econometrics.com/. The toolbox is provided as a zipped ﬁle which can be unzipped to the MATLAB toolbox directory on your PC ( C:\ProgramFiles\MATLAB704\toolbox or my PC - something similar on yours). This should create a subdirectory econometrics in this toolbox directory This econometrics directory will contain a large number of subdirectories containing the various econometric functions. When you next start MATLAB you can access the functions by adding to the path that MATLAB uses to search for functions. You can do the when you next start MATLAB by —File—Set Path— selecting the Add with sub-folders button and navigating to and selecting the econometrics folder. If you select save after entering the directory the functions will be available each time you start MATLAB. If you have the required permissions you can also access the toolbox from the IIS server. The toolbox provides full source code for each function. Thus if no function provides the capabilities that you require it may be possible the amend the function and add the required functionality. If you do such work you should consider submitting you program for inclusion in a future version of the toolbox. By collaborating in this way you are helping to ensure the future of the project The programs in the toolbox are examples of good programming practice and have good comments. If you are starting some serious programming in MATLAB you could learn a lot about programming by reading these programs. Sample run from the LeSage toolbox To illustrate the use of the LeSage toolbox I set out below the output of the demonstration program demo reg.m. This program illustrates many of the various univariate estimation procedures available in the toolbox. % PURPOSE: demo using most all regression functions % % ols,hwhite,nwest,ridge,theil,tsls,logit,probit,tobit,robust %--------------------------------------------------% USAGE: demo_all %--------------------------------------------------clear all; rand(’seed’,10); n = 100; k=3;

64

xtmp = randn(n,k-1); tt = 1:n; ttp = tt’; e = randn(n,1).*ttp; % heteroscedastic error term %e = randn(n,1); b = ones(k,1); iota = ones(n,1); x = [iota xtmp]; % generate y-data y = x*b + e; vnames=strvcat(’yvar’,’iota’,’x1’,’x2’); % * * * * * * * reso = ols(y,x); prt(reso,vnames); % * * * * * * * demo hwhite regression demo ols regression % homoscedastic error term

res = hwhite(y,x); prt(res,vnames); % * * * * * * * demo nwest regression

nlag=2; res = nwest(y,x,nlag); prt(res,vnames); % * * * * * * * demo ridge regresson

rres = ridge(y,x); prt(rres,vnames); % * * * * * * * n = 24; y = zeros(n,1); y(1:14,1) = ones(14,1); % (data from Spector and Mazzeo, 1980) xdata = [21 24 25 26 28 31 33 34 35 37 43 49 ... 51 55 25 29 43 44 46 46 51 55 56 58]; iota = ones(n,1); x = [iota xdata’]; demo logit regression

65

vnames=strvcat(’days’,’iota’,’response’); res = logit(y,x); prt(res,vnames); % * * * * * * * n = 32; k=4; demo probit regression

y = zeros(n,1); % grade variable y(5,1) = 1; y(10,1) = 1; y(14,1) = 1; y(20,1) = 1; y(22,1) = 1; y(25,1) = 1; y(25:27,1) = ones(3,1); y(29,1) = 1; y(30,1) = 1; y(32,1) = 1; x = zeros(n,k); x(1:n,1) = ones(n,1); % intercept x(19:32,2) = ones(n-18,1); % psi variable tuce = [20 22 24 12 21 17 17 21 25 29 20 23 23 25 26 19 ... 25 19 23 25 22 28 14 26 24 27 17 24 21 23 21 19]; x(1:n,3) = tuce’; gpa = [2.66 2.89 3.28 2.92 4.00 2.86 2.76 2.87 3.03 3.92 ... 2.63 3.32 3.57 3.26 3.53 2.74 2.75 2.83 3.12 3.16 ... 2.06 3.62 2.89 3.51 3.54 2.83 3.39 2.67 3.65 4.00 ... 3.10 2.39]; x(1:n,4) = gpa’; vnames=strvcat(’grade’,’iota’,’psi’,’tuce’,’gpa’); resp = probit(y,x); prt(resp,vnames); % results reported in Green (1997, chapter 19) % b = [-7.452, 1.426, 0.052, 1.626 ]

66

% * * * * * * * demo theil-goldberger regression % generate a data set nobs = 100; nvar = 5; beta = ones(nvar,1); beta(1,1) = -2.0; xmat = randn(nobs,nvar-1); x = [ones(nobs,1) xmat]; evec = randn(nobs,1); y = x*beta + evec*10.0; Vnames = strvcat(’y’,’const’,’x1’,’x2’,’x3’,’x4’); % set up prior rvec = [-1.0 1.0 2.0 2.0 1.0]; rmat = eye(nvar); bv = 10000.0; % umat1 = loose prior umat1 = eye(nvar)*bv; % initialize prior variance as diffuse for i=1:nvar; umat1(i,i) = 1.0; end; lres = theil(y,x,rvec,rmat,umat1); prt(lres,Vnames); % * * * * * * * demo two-stage least-squares regression

% prior means for the coefficients

% overwrite diffuse priors with informative prior

nobs = 200;

67

x1 = randn(nobs,1); x2 = randn(nobs,1); b1 = 1.0; b2 = 1.0; iota = ones(nobs,1); y1 = zeros(nobs,1); y2 = zeros(nobs,1); evec = randn(nobs,1); % create simultaneously determined variables y1,y2 for i=1:nobs; y1(i,1) = iota(i,1)*1.0 + x1(i,1)*b1 + evec(i,1); y2(i,1) = iota(i,1)*1.0 + y1(i,1)*1.0 + x2(i,1)*b2 + evec(i,1); end; vname2 = [’y2-eqn ’, ’y1 var ’, ’constant’, ’x2 var ’];

% use all exogenous in the system as instruments xall = [iota x1 x2]; % do tsls regression result2 = tsls(y2,y1,[iota x2],xall); prt(result2,vname2);

% * * * * * * * demo robust regression % generate data with 2 outliers nobs = 100; nvar = 3; vnames = strvcat(’y-variable’,’constant’,’x1’,’x2’); x = randn(nobs,nvar); x(:,1) = ones(nobs,1); beta = ones(nvar,1);

68

evec = randn(nobs,1); y = x*beta + evec; % put in 2 outliers y(75,1) = 10.0; y(90,1) = -10.0; % get weighting parameter from OLS % (of course you’re free to do differently) reso = ols(y,x); sige = reso.sige; % set up storage for bhat results bsave = zeros(nvar,5); bsave(:,1) = ones(nvar,1); % loop over all methods producing estimates for i=1:4; wfunc = i; wparm = 2*sige; % set weight to 2 sigma res = robust(y,x,wfunc,wparm); bsave(:,i+1) = res.beta; end; % column and row-names for mprint function in.cnames = strvcat(’Truth’,’Huber t’,’Ramsay’,’Andrews’,’Tukey’); in.rnames = strvcat(’Parameter’,’constant’,’b1’,’b2’); fprintf(1,’Comparison of alternative robust estimators \n’); mprint(bsave,in); res = robust(y,x,4,2); prt(res,vnames); % * * * * * * * demo regresson with t-distributed errors

res = olst(y,x); prt(res,vnames);

69

% * * * * * * * res = lad(y,x); prt(res,vnames);

demo lad regression

% * * * * * * * demo tobit regression n=100; k=5; x = randn(n,k); x(:,1) = ones(n,1); beta = ones(k,1)*0.5; y = x*beta + randn(n,1); % now censor the data for i=1:n if y(i,1) < 0 y(i,1) = 0.0; end; end; resp = tobit(y,x); vnames = [’y ’,

’iota ’, ’x1var ’, ’x2var ’, ’x3var ’, ’x4var ’]; prt(resp,vnames);

% * * * * * * * demo thsls regression clear all; nobs = 100; neqs = 3; x1 = randn(nobs,1); x2 = randn(nobs,1); x3 = randn(nobs,1); b1 = 1.0; b2 = 1.0;

70

b3 = 1.0; iota = ones(nobs,1); y1 = zeros(nobs,1); y2 = zeros(nobs,1); y3 = zeros(nobs,1); evec = randn(nobs,3); evec(:,2) = evec(:,3) + randn(nobs,1); % create cross-eqs corr % create simultaneously determined variables y1,y2 for i=1:nobs; y1(i,1) = iota(i,1)*10.0 + x1(i,1)*b1 + evec(i,1); y2(i,1) = iota(i,1)*10.0 + y1(i,1)*1.0 + x2(i,1)*b2 + evec(i,2); y3(i,1) = iota(i,1)*10.0 + y2(i,1)*1.0 + x2(i,1)*b2 + x3(i,1)*b3 + evec(i,3); end;

vname1 = [’y1-LHS ’, ’constant’, ’x1 var vname2 = [’y2-LHS ’y1 var ’]; ’, ’,

’constant’, ’x2 var ’]; vname3 = [’y3-LHS ’y2 var ’, ’,

’constant’, ’x2 var ’, ’x3 var ’];

% set up a structure for y containing y’s for each eqn y(1).eq = y1; y(2).eq = y2; y(3).eq = y3; % set up a structure for Y (RHS endogenous) for each eqn Y(1).eq = []; Y(2).eq = [y1]; Y(3).eq = [y2];

71

% set up a structure fo X (exogenous) in each eqn X(1).eq = [iota x1]; X(2).eq = [iota x2]; X(3).eq = [iota x2 x3]; % do thsls regression result = thsls(neqs,y,Y,X); vname = [vname1 vname2 vname3]; prt(result,vname); % * * * * * * * demo olsc, olsar1 regression % generate a model with 1st order serial correlation n = 200; k = 3; tt = 1:n; evec = randn(n,1); xmat = randn(n,k); xmat(:,1) = ones(n,1); beta = ones(k,1); beta(1,1) = 10.0; % constant term y = zeros(n,1); u = zeros(n,1); for i=2:n; u(i,1) = 0.4*u(i-1,1) + evec(i,1); y(i,1) = xmat(i,:)*beta + u(i,1); end; % truncate 1st 100 observations for startup yt = y(101:n,1); xt = xmat(101:n,:); n = n-100; % reset n to reflect truncation Vnames = [’y ’,

’cterm’,

72

’x2 ’x3

’, ’];

% do Cochrane-Orcutt ar1 regression result = olsc(yt,xt); prt(result,Vnames); % do maximum likelihood ar1 regression result2 = olsar1(yt,xt); prt(result2,Vnames);

% * * * * * * * demo switch_em, hmarkov_em regressions clear all; % generate data from switching regression model nobs = 100; n1 = 3; n2 = 3; n3 = 3; b1 = ones(n1,1); b2 = ones(n2,1)*5; b3 = ones(n3,1); sig1 = 1; sig2 = 1; randn(’seed’,201010); x1 = randn(nobs,n1); x2 = randn(nobs,n2); x3 = randn(nobs,n3); ytruth = zeros(nobs,1); for i=1:nobs; if x3(i,:)*b3 0 d= 0 if y ≤ 0. The log-likelihood can then be written
N

i=1

1 1 1 di (− ln 2π − ln σ 2 − 2 (yi − xi β)2 ) + (1 − di ) ln 1 − Φ 2 2 2σ

xβ σ

1. Matlab program to calculate tobit log-likelihood This is an amended version of a sample program included with the Le Sage econometrics package function like = to_liked(b,y,x); % PURPOSE: evaluate tobit log-likelihood % to demonstrate optimization routines %----------------------------------------------------% USAGE: % where: % like = to_liked(b,y,x) b = parameter vector (k x 1) y = dependent variable vector (n x 1)

% x = explanatory variables matrix (n x m) %----------------------------------------------------% NOTE: this function returns a scalar equal to the % negative of the log likelihood function % % or a scalar sum of the vector depending on the value of the flag argument

% k ~= m because we may have additional parameters % in addition to the m bhat’s (e.g. sigma) %----------------------------------------------------% error check if nargin ~= 3,error(’wrong # of arguments to to_like1’); end; [m1 m2] = size(b); if m1 == 1 b = b’; end;

83

h = .000001; [m junk] = size(b); beta = b(1:m-1); sigma = max([b(m) h]);

% avoid sigma = 0 % pull out bhat % pull out sigma

xb = x*beta; llf1 = -0.5*log(2*pi) - 0.5*log(sigma^2) -((y-xb).^2)./(2*sigma^2); %amended xbs = xb./sigma; cdf = .5*(1+erf(xbs./sqrt(2))); %amended llf2 = log(h+(1-cdf)); llf = (y > 0).*llf1 + (y

Similar Documents

Forecasting Models

Human Step Detection and Counting with a Microphone (June 2015)

Matlab

Matlab

Matlab

Matlab

Matlab

Matlab

Matlab

Matlab Ode

Matlab

Matlab & Ss

Pelajaran Matlab

Ocr Matlab

Matlab Sol