ECLR - User contributions [en]

Temp

2017-10-06T12:07:27Z

AS: AS moved page Monte-Carlo to Temp: AH request

= Monte Carlo Basics =

Quite often we would like to know which of the proposed estimators is “better” and under what conditions. The exact answer (for any sample size) to these questions is known only for a quite limited class of estimators, distributions of error terms and regression functional forms. We can address these questions using various Central Limit Theorems. However, the applicability of our results relies on a “sufficiently large” number of observations for asymptotics to kick in. The question “When ''sufficiently large'' is ''large enough''?” is a very interesting one, but it lies well beyond the scope of this discussion. In small samples asymptotics does not work. Unfortunately, this is exactly the area where all real-world applications belong to.

Monte Carlo method can somewhat address this issue. However, it has its own problems. Assume that we would like to know how the conditional quasi-maximum likelihood estimator (estimator A) performs compared to the unconditional quasi-maximum likelihood estimator (estimator B) for AR(1) process with AR coefficient <math>\alpha=.5</math>, sample size <math>T=50</math> and t-distributed error terms with 10 degrees of freedom. On one hand, using MC method, we can reliably test the hypotheses that, for example, estimator A has smaller bias and/or smaller variance than estimator B. However, this answer, strictly speaking, can be applied '''only''' for the case of <math>\alpha=.5</math>, <math>T=50</math>, and t-disributed error term with 10 degrees of freedom. Once the sample size, autoregressive coefficient, degrees of freedom or distribution is changed, results might change significantly as well.

The idea behind the Monte Carlo Method is very simple. Assume for a moment that you can draw somehow observations from an (un)known distribution of interest. Then, by drawing sufficiently large number of observations from this distribution, you can approximate with an arbitrary level of precision:

# finite population moments of any functions of random variables
# CDF and/or pdf and thus quantiles of interest

As an example, for the case when the distribution is known and random number generator is readily available, we can compute

<math>E(|X|),E(|X|^3) \text{ for } X\sim N(0,1)</math>

Formally, taking integrals we’ll get:

<math>\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty} |x|e^{-x^2/2}dx=\frac{2}{\sqrt{2\pi}}\int_{0}^{\infty} e^{-x^2/2}d(x^2/2)=\sqrt{\frac{2}{\pi}}\approx 0.7979</math>

<math>\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty} |x|^3e^{-x^2/2}dx=\frac{2}{\sqrt{2\pi}}\int_{0}^{\infty} x^2e^{-x^2/2}d(x^2/2)=2\sqrt{\frac{2}{\pi}}\approx 1.5958</math>

Alternatively, the Monte-Carlo algorithm to estimate <math>E(|X|)</math> and <math>E(|X|^3)</math> is:

# Draw a large number (say, 10000000) of observations from a normal distribution <math>N(0,1)</math> and compute <math>|x|</math> and <math>|x|^3</math> for each of them.
# Take averages and report results

<source>rng(1);x=randn(10000000,1);
absx=abs(x);
absx3=absx.^3;
disp([mean(absx) mean(absx3)])
0.7982 1.5973</source>
Which is quite close to the analytically evaluated ones. The more observations you take, the closer your results to the analytical ones will be.

However, in reality, the random number generator for the distribution of interest is not readily available. Quite often you have to construct it first.

= OLS properties of AR(1) process estimator =

Consider a very simple example: (un)biasedness of OLS estimator while applying it to AR(1) process

<math>y_t=\alpha y_{t-1}+e_t,\ e_t\sim N(0,\sigma^2)</math>

Recall that the OLS estimator for regression

<math>y=X\beta+e, e\sim (0,\sigma^2_e I)</math>

where <math>I</math> is an identity matrix, is <math>\hat \beta=(X'X)^{-1}X'y=\beta+(X'X)^{-1}X'e</math> Unbiasedness of OLS estimator, assuming <math>(X'X)^{-1}</math> exists, is satisfied as long as <math>E[(X'X)^{-1}X'e]=0</math>. Sufficient condition for this to hold is <math>E(e|X)=0</math>. While this is quite often the case for cross-sectional applications, this is definitely not the case for Time Series. The condition <math>E(e|X)=0</math> has to be understood in the following way: ''any'' error term <math>e_s</math> in this regression has to be independent from ''any'' <math>x_t</math>. For AR processes, for obvious reasons, this is not the case. For example, <math>e_t</math> is not independent from <math>x_{t+1}</math> as <math>x_{t+1}</math> is <math>y_t</math>. <math>Cov(x_{t+1},e_t)=\sigma_e^2</math>. Moreover, <math>Cov(x_{t+s},e_t)=\alpha^s\sigma_e^2</math> and, thus the larger <math>|\alpha|</math>, the larger, presumably, has to be the bias. On the other hand, for consistency, it is sufficient that

<math>\label{consistency}
(X'X/T)^{-1} \Rightarrow Q,\ X'e/T \Rightarrow 0</math>

• where <math>Q</math> is a full-rank finite matrix. Since these conditions are satisfied for covariance-stationary processes with gaussian increments, we expect that the bias will become smaller and smaller as the sample size increases.

Let’s sketch the algorithm:

# Generate <math>y_t</math> for <math>t=1,\ldots,T</math>, where <math>T=20</math> and <math>\alpha=0</math> using the following data generating process: <math>y_t=\alpha y_{t-1}+e_t</math>
# Estimate and save <math>\alpha</math> using OLS (Please note, these steps correspond to the first step from the previous algorithm “draw a large number of observations from a distribution of interest)
# Repeat pp 1-2 1000 times, compute <math>\bar{\hat{\alpha}}</math> and bias
# Repeat pp 1-3 for <math>\alpha=0.1, 0.2, ... , 0.9, 0.99</math>
# Repeat pp 1-4 for <math>T=50,100,200,1000</math>
# Plot results for different sample sizes

Before implementing the algorithm, we have to take care of one additional issue, i.e. sample simulation (p. 1). If we know <math>y_0</math>, we can generate <math>y_1</math> drawing an observation from a normal distribution (using normal random number generator <code>randn</code>) and proceed to <math>y_T</math>.

<math>\begin{aligned}
y_1&=\alpha y_{0}+e_1\\
y_2&=\alpha y_{1}+e_2\\
\ldots\\
y_T&=\alpha y_{T-1}+e_T\end{aligned}</math>

• Unfortunately, <math>y_0</math> is unknown. ''Please note: The whole discussion regarding <math>y_0</math> below applies for processes with finite memory only. Compare with Random Walk.''

One way to handle this problem is ''to pretend that <math>y_0</math> is known'', i.e. to assume that <math>y_0=E(y)=0</math>. Then, <math>y_1=e_1</math>, <math>y_2=\alpha e_1+e_2</math>, etc. However, this approach has an additional issue. For AR(1) process <math>Var(y_t)=\frac{\sigma_e^2}{1-\alpha^2}</math> for all <math>t</math>. If we assume that <math>y_0=0</math>, however, <math>Var(y_1)=Var(e_1)=\sigma_e^2</math>, <math>Var(y_2)=\sigma_e^2(1+\alpha^2)</math>, and for arbitrary <math>t>0</math>: <math>Var(y_t)=\sigma_e^2\sum_{i=0}^{t-1}\alpha^{2(t-1)}=\sigma^2_e\frac{1-\alpha^{2t}}{1-\alpha^2}</math>. For a reasonably large <math>t</math>: <math>V(y_t)\approx\frac{\sigma^2_e}{1-\alpha^2}</math>. Thus, if we decide to implement this algorithm, we have to start recording <math>y</math> not from the first observation, but from some non-zero <math>T_0</math> s.t. <math>\alpha^{2T_0}\approx 0</math>. Obviously, <math>T_0</math> depends on <math>\alpha</math>. The larger <math>\alpha</math>, the larger <math>T_0</math> should be. The thrown away part is sometimes called a “burn-in” period.

Alternatively, we can take into account the fact that <math>y_0</math> is not known and generate <math>y_1</math> directly with <math>E(y_1)=0</math> and <math>Var(y_1)=\frac{\sigma^2_e}{1-\alpha^2}</math>, i.e. <math>y_1=e_1/\sqrt{1-\alpha^2}</math>

For stationary processes the first approach does not require any analytical derivations, i.e. any knowledge of unconditional and various conditional population moments of the process. The second one runs faster. However, the second method has its own issues. It is easy to implement it for AR(1) process. For ARMA(p,q) the amount of effort increases significantly. For stochastic volatility models it might become unbearable.

I will focus on the first version. <math>T_0=200</math> is sufficiently large for almost all values of <math>\alpha</math> considered. The p.1 of the algorithm becomes:

# Generate <math>y_t</math> for <math>t=1,\ldots,T</math>, where <math>T=20</math> and <math>\alpha=0</math>
## Draw <math>T+T_0</math> observations from <math>N(0,1)</math> with <math>T_0=200</math>, save them in a vector <math>e</math> and assume <math>z_0=0</math>
## Construct <math>z_t</math> for <math>t=1</math>: <math>z_t=\alpha z_{t-1}+e_t</math>
## Repeat 1.b) for <math>t=2,3,\ldots,T+T_0</math>
## Save the last T observations as vector <math>y</math>

To make the code easier to understand, we will write a set of functions for these steps. (Please see the section about functions). Also, we will use <code>for .. end</code> loop here. Please see the section about loops.

Lets implement the first step, i.e. generation of a sample of size <math>T</math>. This function has to have three parameters: alpha, <math>T_0</math>, <math>T</math>. We call it dgp, i.e. data generating process.

<source>function y=dgp(alpha,T0,T)
Tt=T+T0;
e=randn(Tt,1); %1 a)
z=zeros(Tt,1);
z(1)=e(1); %since z_0=0
for t=2:Tt
z(t)=alpha*z(t-1)+e(t);
end
y=z(end-T+1:end); %dropping first T0 observations
end</source>
The next step is to construct an estimate for <math>\alpha</math>:

<source>function alpha=myar1(y)
X=y(1:end-1);%specifying X
y=y(2:end);%specifying y
alpha=inv(X'*X)*X'*y; %Estimating AR(1) parameter
end</source>
The next step is to repeat these steps 1000 times, compute mean, bias and standard errors of <math>\hat\alpha</math>. For the sake of generality, this function will depend on <math>T_0</math>, <math>T</math>, <math>\alpha</math>, <math>N</math>, where <math>N</math> - number of replication, in our case 1000.

<source>function [meanalpha biasalpha stdalpha]=MC(alpha,T0,T,N)
for i=N:-1:1
%This way we don't need to initialize alphahat vector prior to the cycle
y=dgp(alpha,T0,T);
alphahat(i,1)=myar1(y);
end
meanalpha=mean(alphahat);
biasalpha=mean(alphahat-alpha);
stdalpha=std(alphahat);
end</source>
The next step is to repeat our simulations for different sample sizes and different alphas. The easiest way is to create two vectors of values, one for different sample sizes, one for different alphas, and loop over all possible combinations of them:

<source>alphavec=[0:.1:.9 .99];
Tvec=[20 50 100 200 1000];
T0=200; %burn-in period
N=1000; %number of iterations
rng(1); %for results replicability
for Tind=1:length(Tvec)
T=Tvec(Tind); %picking the right T
for alphaind=1:length(alphavec)
alpha=alphavec(alphaind); %picking the right alpha
[meanalpha(alphaind,Tind) biasalpha(alphaind,Tind) stdalpha(alphaind,Tind)]=MC(alpha,T0,T,N);
end
end</source>
And the last, but not the least, step is to plot our results:

<source>plot(alphavec,biasalpha)
title('Bias as a function of sample size and $\alpha$','interpreter','latex')
xlabel('$\alpha$','interpreter','latex')
ylabel('Bias of $\hat\alpha$','interpreter','latex')
legend({'T=20','T=50','T=100', 'T=200', 'T=1000'})
figure
plot(alphavec,stdalpha)
title('Standard error of $\hat\alpha$ as a function of sample size and $\alpha$','interpreter','latex')
xlabel('$\alpha$','interpreter','latex')
ylabel('Std of $\hat\alpha$','interpreter','latex')
legend({'T=20','T=50','T=100', 'T=200', 'T=1000'})</source>
The next two figures are generated using the code above. To create a smoother picture, we ran simulations for <math>N=10000</math>, <math>\alpha=0,0.05,0.1,0,15,\ldots,0.95,0.99</math>

[[Image:bias.png|image|thumb|center|upright=3.0|alt=Bias of AR(1) Estimation.|Dependence of Bias on $\alpha$.]]

[[Image:std.png|image|thumb|center|upright=4.0|alt=St. deviation of AR(1) Estimation.|Dependence of std on $\alpha$.]]

Some of the findings coincide with our predictions, some of them do not. In particular,

# Bias of <math>\alpha=0</math> is very close to 0
# Bias of <math>\alpha</math> increases (in absolute value) up to some point
# Standard errors decrease (Figure 2) as sample size increases.

One outcome that contradicts our “intuition” is the fact that bias in absolute value ''decreases'' after some point. This point depends on the sample size. The intuition for this behavior lies in the notion of non-stationarity. As <math>\alpha</math> further and further increases, stationary series start to behave more and more like non-stationary ones. For non-stationary time series <math>X'X/T^2=\sum y_t^2 /T^2 \Rightarrow C(B(t))</math>, while <math>X'e/T = \sum y_{t-1}e_t/T \Rightarrow C'(B(t))</math>, where <math>C(B(t)),C'(B(t))</math> – finite functionals of Brownian motions. Thus, for <math>|\rho|\rightarrow 1</math>, <math>(X'X)^{-1}</math> will dominate <math>X'y</math> and the bias will decrease. Moreover, for <math>\rho=1</math>, <math>\hat\alpha</math> is a consistent estimator of <math>\alpha</math> even in the case of <math>Cov(x_t,e_t)\ne 0</math>.

= On a notion of Conditional Distribution =

Repeat simulation using a slightly modified version of the main code, i.e.

<source>alphavec=[0:.1:.9 .99];
Tvec=[20 50 100 200 1000];
T0=200; %burn-in period
N=1000; %number of iterations
%rng(1); %for results replicability
for Tind=1:length(Tvec)
T=Tvec(Tind); %picking the right T
for alphaind=1:length(alphavec)
rng(1); %this is the new line
alpha=alphavec(alphaind); %picking the right alpha
[meanalpha(alphaind,Tind) biasalpha(alphaind,Tind) stdalpha(alphaind,Tind)]=MC(alpha,T0,T,N);
end
end</source>
Compare the smoothness of the graphs. Explain the differences.

Monte-Carlo

2017-10-06T12:07:27Z

AS: AS moved page Monte-Carlo to Temp: AH request

#REDIRECT [[Temp]]

MATLAB

2017-10-06T11:14:43Z

AS:

== <div id="Essential"></div>The Essential MATLAB Programming Techniques ==

In this section we will introduce a number of basic and intermediate programming techniques. Whatever language you program in you will encounter these techniques, although the details will, of course, vary. We recommend that you ensure that you are familiar with these before you progress to [[#SpecEcmtrTopics| Special Econometric Topics ]].

=== Basic Programming ===

{| class="wikitable"
|-
! scope="col"| Basics and Matrices
! scope="col"| Loading Data and Date Formats
! scope="col"| Program Flow and Logicals
! scope="col"| Functions
! scope="col"| Saving Data and Screen Output
|-
| [[Discussion]] [http://www.youtube.com/watch?v=av5MgVpybT0&feature=youtu.be&hd=1 Example Clip]
| [[LoadingData|Discussion]] [http://youtu.be/jyb68zGM2ik?hd=1 ExampleClip]
| [[Program Flow and Logicals|Discussion]]
| [[Function|Discussion]] [[FctExampleCode|Example Code]] [[media:OLSexample.xls|OLSexample.xls]] [http://youtu.be/FPw9DH8pfiU?hd=1 Example Clip]
| [[SavingData|Discussion]]
|}

After having gone through these basic techniques you may want to test your newly acquired skills with the following examples.

{| class="wikitable"
|-
! scope="col"| Example 1
! scope="col"| Example 2
|-
| [[Example 1]]
| [[Example 2|Example2a]] [[Example 2b|Example2b]]
|}

=== Intermediate Programming ===

{| class="wikitable"
|-
! scope="col"| Statistical Functions
! scope="col"| Arrays and Structures
! scope="col"| Debugging
! scope="col"| Graphing Data
! scope="col"| Function Handlers Anonymous Functions
|-
| [[StatFunct|Discussion]]
| [[ArrayStructures|Discussion]]
| coming soon
| [[Graphing|Discussion]]
| [[Anonym|Discussion]]
|}

=== Advanced Programming ===

Sorry, but this cannot be taught! It will come with experience. Find someone who has experience in MATLAB programming and let him or her look over your code.

== Nonlinear Optimisation ==

The optimal parameters in a linear econometric model (assuming certain assumptions) can be found analytically. We call them the Ordinary Least Squares (OLS) estimates and they are easily calculated with a certain formula (see the [[FctExampleCode#OLSestm|OLSest.m]] function). When econometric models do not have such an analytical solution, an alternative parameter estimation strategy is required. In essence it is a clever "trial and error" strategy. This is often called nonlinear optimisation.

Nonlinear optimisation is a very important, but also a very tricky area of econometric computing. It certainly helps to understand some of the underlying theory and therefore we have below separate sections on the theory and implementation.

{| class="wikitable"
|-
! scope="col"| Theory
! scope="col"| Implementation
! scope="col"| Constrained Optimisation
|-
| [[NonlinOptTheory| Discussion]]
| [[NonlinOptImp| Discussion]]
| [[ConNonlinOptImp| Discussion]]
|}

== <div id="SpecEcmtrTopics"></div> Special Econometric Topics ==

Topics in this Section will assume that you have mastered all the techniques covered in the [[#Essential| Essential Programming Section ]]

{| class="wikitable"
|-
! scope="col"| Robust standard errors
! scope="col"| Univariate Time Series
! scope="col"| Unit Root and Stationarity Testing
! scope="col"| Forecasting
|-
| [[RobInf|Discussion]] [[ExampleCodeOLShac|Example Code]]
| [[UniTS|Discussion]] [[media:FXrateUSEU.xls|FXrateUSEU.xls]] [[media:USGDP.xls|USGDP.xls]]
| coming soon
| [[Forecasting|Discussion]]
|}

{| class="wikitable"
|-
! scope="col"| Maximum Likelihood
! scope="col"| Generalized Methods of Moments
! scope="col"| Instrumental Variables
! scope="col"| Bayesian Estimation
|-
| [[MaxLik|Discussion]] [[MaxLikCode|Example Code]]
| coming soon
| [[IV|Discussion]] [[ExampleCodeIV|Example Code]]
| [[Bayes|Discussion]]
|}

{| class="wikitable"
|-
! scope="col"| Monte-Carlo/ Simulation Techniques
! scope="col"| Binary Response Models
! scope="col"| Handling High Frequency Data
|-
| coming soon
| coming soon
| coming soon

|}

== Other useful MATLAB resources ==

=== The MATLAB Software ===

The software is available on University of Manchester Computer Labs. If you make regular use of MATLAB you should consider purchasing your own software. The Student Version of MATLAB is available, for instance, from [http://www.amazon.co.uk/MATLAB-Simulink-Student-Version-R2014a/dp/0989614026/ref=sr_1_1?s=software&ie=UTF8&qid=1411983990&sr=1-1&keywords=matlab+2014 Amazon] for £66. This is a real bargain, considering that the equivalent non-discounted package would come in at about £4,000.

=== Freely available toolboxes ===

The following toolboxes are freely available and contain extremely useful procedures

* Spatial Econometrics by James P. LeSage [http://www.spatial-econometrics.com/]. This toolbox contains a wide variety of useful econometrics functions. It also contains an excellent documentation. In addition to quite general econometric functions you will, as the name suggests, find a huge list of functions relevant if you are working with spatial data.
* <div id="MFEtoolbox"></div>Oxford MFE toolbox by Kevin Sheppard [https://bitbucket.org/kevinsheppard/mfe_toolbox]. Use the download link in the box on the right that starts with "Owner: Kevin Sheppard". This toolbox contains many useful functions for uni- and multivariate volatility models.

You need to copy these toolboxes into your MATLAB toolbox folder and add the respective path to the MATLAB list of folders it searches for functions. (In the main menu select FILE and then SET PATH where you can add the folders you added.) If you work on a computer for which you have no administrator rights, this strategy may not work. This [http://youtu.be/_32OqcW9WoY?hd=1 Example Clip] demonstrates what to do in that case. It is just a matter of adding one line into your code! Piece of cake.

=== Literature and other learning resources ===
* [http://www.kevinsheppard.com/wiki/MFE_Toolbox: Kevin Sheppard's MATLAB introduction].
* Martin V., Hurn S. and Harris D. (2012) Econometric Modelling with Time Series: Specification, Estimation and Testing (Themes in Modern Econometrics).[http://www.amazon.co.uk/Econometric-Modelling-Time-Specification-Econometrics/dp/0521196604/ref=sr_1_1?s=books&ie=UTF8&qid=1345214275&sr=1-1] This book contains an extensive library of relevant MATLAB codes.
* Higham, D.J. and Higham, N.J. (2005) MATLAB Guide, Society for Industrial and Applied Mathematics [http://www.amazon.co.uk/MATLAB-Guide-Desmond-J-Higham/dp/0898715784/ref=sr_1_1?s=books&ie=UTF8&qid=1347377409&sr=1-1]
This website does not cover any theoretical ground and is no substitute for any Econometric Textbook. There is a wide range of very good Econometric Textbooks available. If you are concerned about programming in MATLAB than you are likely to appreciate textbooks that use matrix notation. Here are two very good books that fit that bill:
* Heij C., de Boer P., Franses P.H., Kloek T. and van Dijk H.K (2004) Econometric Methods with Applications in Business and Economics, Oxford University Press, New York.[http://www.amazon.co.uk/Econometric-Methods-Applications-Business-Economics/dp/0199268010/ref=sr_1_1?s=books&ie=UTF8&qid=1354473313&sr=1-1]
* Greene W.H. (2012) Econometric Analysis, Pearson, Harlow.[http://www.amazon.co.uk/Econometric-Analysis-William-H-Greene/dp/0273753568/ref=sr_1_1?ie=UTF8&qid=1354473593&sr=8-1]

Program Flow and Logicals

2013-10-11T13:25:57Z

AS: /* while end loop */

= Preliminaries =

Very often in your life you have to repeat the same operation many times (move your right and left leg sequentially while walking/running) or behave differently depending on external conditions (there is or there isn’t a bus on the bus stop). Quite often these two are combined together. Say, if there is a bus on the bus stop, then you run trying to catch it, otherwise walk or stop and enjoy the usual Manchester weather. The same is true for programming. Quite often you want to repeat the same operation many times, or you want to change the way you process your data depending on some conditions. We start with conditional statements. They execute different pieces of the code depending on whether <source enclose=none>condition</source> is true or false. There are several ways you can formulate it. The shortest

<source>if condition
statement1;
statement2;
...
end
</source>
, executes <source enclose=none>statement1;statement2,...</source> only if <source enclose=none>condition</source> is satisfied. <source enclose=none>condition</source> can be anything that generate non-zero or 0 (True or False), say <source enclose=none>i>0</source>, <source enclose=none>size(y,1)\sim=40</source>, or <source enclose=none>5-i</source>. The last condition is True always but for <source enclose=none>i=5</source>. A slightly longer version

<source> if condition
statement1;
statement2;
...
else
statement1a;
statement2a;
...
end
</source>
runs <source enclose=none>statement1;statement2;...</source>, if <source enclose=none>condition</source> is true and <source enclose=none>statement1a;statement2a;...</source> otherwise. The most general specification is

<source> if condition1
statement1;
statement2;
...
elseif condition2
statement1a;
statement2a;
...
...
...
elseif conditionN
statement1b;
statement2b;
...
else
statement1c;
statement2c;
...
end</source>
In this case, however, you have to ensure that <source enclose=none>condition1, condition2, …, conditionN</source> are mutually disjoint. As an example, you might think about different actions depending on your final grade. Condition1: <source enclose=none>grade<30</source>; Condition 2: <source enclose=none>(grade>=30)&&(grade<40)</source>; Condition 3: <source enclose=none>(grade>=40)&&(grade<50)</source>; etc.

MATLAB has two statements that create a loop. First, it is unconditional loop:

<source>for CounterVariable=[range of values]
statement1;
statement2;
...
end</source>
It repeats at most as many times as many elements it has in the <source enclose=none>[range of values]</source>. If the range of values is empty, this loop does not run. Say, if you define a range <source enclose=none>10:1</source>, MATLAB creates an empty range. Thus, this loop will not be executed. If you define a range <source enclose=none>1:3:10</source>, MATLAB creates a range of four values <source enclose=none>[1 4 7 10]</source>, and the loop runs four times. During the first iteration, <source enclose=none>CounterVariable=1</source>, during the second <source enclose=none>CounterVariable=4</source>, etc. After the end of the loop <source enclose=none>CounterVariable=10</source>. Please note, it is very unwise to modify the counter inside the loop. All modifications will disappear after the next iteration. Please also note, that the values in the range could be anything, including filenames or matrices from a cell vector. There are two commands that can modify the execution of the loop. <source enclose=none>continue</source> breaks the current ''iteration'' of the loop. Once it is executed, the loop continues skipping current iteration. <source enclose=none>break</source> stops the execution of the loop and your program continues after this point. These commands are used inside <source enclose=none>if</source> statements. For example, <source enclose=none>if CounterVariable == 10 continue;end</source> skips the loop iteration for <source enclose=none>CounterVariable = 10</source>.

Second, it is the conditional loop

<source>while condition
statement1;
statement2;
...
end</source>
This version of the loop executes statements as long as <source enclose=none>condition</source> is true. If <source enclose=none>condition</source> is always true, your loop runs forever.

== <source enclose=none>for ... end</source> loop ==

A standard application for <source enclose=none>for ... end</source> loop is the reconstruction of AR(p) series once AR(p) coefficients and the vector of error terms is known. <math>y_t=\phi_0+\sum_{i=1}^p \phi_i y_{t-i}+e_t.</math> For simplicity, we assume that <math>p=1</math>. Also, to be able to compute <math>y_1</math>, we need to provide <math>y_0</math>. Since we don’t know <math>y_0</math>,the best guess for <math>y_0</math> is <math>E(y_0)</math>. For stationary AR(1) process, that is for the case <math>|\phi_1|<1</math>, <math>E(y_0)=\phi_0/(1-\phi_1)</math>. Thus, knowing <math>y_0</math> and <math>e_t</math> for <math>t=1,\ldots,T</math>, we can reconstruct <math>y_t,\ t=1\ldots,T</math>:

<math>\begin{aligned}
y_1=&\phi_0+\phi_1 y_0+e_1\\
y_2=&\phi_0+\phi_1 y_1+e_2\\
&\ldots\\
y_t=&\phi_0+\phi_1 y_{t-1}+e_t\\
&\ldots\\
y_T=&\phi_0+\phi_1 y_{T-1}+e_T\end{aligned}</math>

Definitely, if you are patient enough and <math>T</math> is not very large, you can create your m file with <math>T</math> lines in it. However, once <math>T</math> is unknown, this approach would not work. Fortunately, there is a better alternative for this type of operations. All these computations can be summarized using the following algorithm, assuming vector <source enclose=none>e</source> is already created:

# Find the length of a vector of error terms <source enclose=none>e</source>: <source enclose=none>T=size(e,1)</source>
# Initialize a vector <source enclose=none>y</source> of the same length as vector <source enclose=none>e</source>: <source enclose=none>y=zeros(T,1)</source>
# Compute <source enclose=none>y(1)=phi0+phi1*(phi0/(1-phi1))+e(1)</source>. Please remember, we assume that <math>y_0=E(y)=\phi_0/(1-\phi_1)</math>
# Compute <source enclose=none>y(i)=phi0+phi1*y(i-1)+e(i)</source> for <math>i=2</math>
# Repeat line 4 for <math>i=3,...,T</math>

When vector <source enclose=none>e</source> is known in advance (for the sake of certainty, let <source enclose=none>e=randn(1000,1)</source>), the MATLAB code is
<source>
T=size(e,1);
y=zeros(T,1);
y0=phi0/(1-phi1);
y(1)=phi0+phi1*y0+e(1);
for i=2:T
y(i)=phi0+phi1*y(i-1)+e(i);
end</source>
However, if <source enclose=none>phi1=1</source>, <math>E(y_t)</math> is not constant, then the formula we use in the code does not work and will create either a series <source enclose=none>y</source> of <math>\pm\infty</math>, if <source enclose=none>phi0 \ne 0</source> or a series of not a numbers <source enclose=none>NaN</source>, if <source enclose=none>phi0=0</source><ref>There are two special numerical values in MATLAB. One is infinity <source enclose=none>Inf</source>, and another is not a number <source enclose=none>NaN</source>. A value of a variable becomes <source enclose=none>Inf</source> if the number is too big in absolute value (<math>\approx \pm 2e308</math>). Also, infinity is generated once you have expressions like <math>x/0</math>, where <math>x\ne0</math>. After that, infinity can only change a sign or become not a number. Not a number appears when there is an uncertainty of a kind of <math>0/0</math>, <math>\infty-\infty</math> and such. Any algebraic operations with <source enclose=none>NaN</source> result <source enclose=none>NaN</source>
</ref>.

== <source enclose=none>if else end</source> or <source enclose=none>if end</source> ==

To avoid these inconveniences, we have to consider separately two cases:

# AR(1) process is stationary, i.e. <math>|\phi_1|<1</math>
# AR(1) process is nonstationary, i.e. <math>|\phi_1|\ge 1</math>

For the latter, we have to acknowledge the fact that <math>E(y_t)=\mu_t</math>, i.e. unconditional expectation is a function of time. In this case we have to set <math>E(y_0)</math> to some value. A standard assumption for non-stationary series is to assume that <math>E(y_0)=0</math>.

The algorithm in this case would look like:

# Find a length of a vector of error terms <source enclose=none>e</source>: <source enclose=none>T=size(e,1)</source>
# Initialize a vector <source enclose=none>y</source> of the same length as vector <source enclose=none>e</source>: <source enclose=none>y=zeros(T,1)</source>
# Check whether <source enclose=none>abs(phi1)<1</source>. If this statement is true, then <source enclose=none>y0=phi0/(1-phi1)</source>. Else, <source enclose=none>y0=0</source>. Please remember, we set <math>y_0=E(y_0)</math>.
# Compute <source enclose=none>y(1)=phi0+phi1*y0+e(1)</source>.
# Compute <source enclose=none>y(i)=phi0+phi1*y(i-1)+e(i)</source> for <math>i=2</math>
# Repeat line 5 for <math>i=3,...,T</math>

Assuming vector <source enclose=none>e</source> is known in advance (for the sake of certainty, let <source enclose=none>e=randn(1000,1)</source>), the MATLAB code is

<source> T=size(e,1);
y=zeros(T,1);
if abs(phi1)<1
y0=phi0/(1-phi1);
else
y0=0;
end
y(1)=phi0+phi1*y0+e(1)
for i=2:T
y(i)=phi0+phi1*y(i-1)+e(i);
end</source>
If you don’t like the word else, you can skip it:

<source> T=size(e,1);
y=zeros(T,1);
y0=0;
if abs(phi1)<1
y0=phi0/(1-phi1);
end
y(1)=phi0+phi1*y0+e(1)
for i=2:T
y(i)=phi0+phi1*y(i-1)+e(i);
end</source>

== <source enclose=none>while end</source> loop ==

An alternative way of running the same code is to use a conditional loop (purely for demonstration purposes). Usually the conditional loop is used when the number of iterations is not known in advance.

# Find the length of a vector of error terms <source enclose=none>e</source>: <source enclose=none>T=size(e,1)</source>
# Initialize a vector <source enclose=none>y</source> of the same length as vector <source enclose=none>e</source>: <source enclose=none>y=zeros(T,1)</source>
# Check whether <source enclose=none>abs(phi1)<1</source>. If this statement is true, then <source enclose=none>y0=phi0/(1-phi1)</source>. Else, <source enclose=none>y0=0</source>. Please remember, we set <math>y_0=E(y_0)</math>.
# Set $i=1$.
# Compute <source enclose=none>y(i)=phi0+phi1*y0+e(i)</source>.
# Increase $i$ by one, i.e. $i=i+1$ (please note, in programming this is a valid statement)
# Check whether <math>i<=T</math>. If yes, proceed. Else go to line 11.
# Compute <source enclose=none>y(i)=phi0+phi1*y(i-1)+e(i)</source>
# Increase i by 1, i.e. <math>i=i+1</math>
# Repeat lines 7 $-$ 9
# ... (continuing the script)
Assuming vector <source enclose=none>e</source> is known in advance (for the sake of certainty, let <source enclose=none>e=randn(1000,1)</source>), the MATLAB code is

<source> T=size(e,1);
y=zeros(T,1);
if abs(phi1)<1
y0=phi0/(1-phi1);
else
y0=0;
end
y(1)=phi0+phi1*y0+e(1)
i=2;
while i<=T
y(i)=phi0+phi1*y(i-1)+e(i);
i=i+1;
end</source>

== Imperfect substitutes of the above ==

MATLAB has two powerful tools that make programmer’s life much easier and utilization of loops/if less frequent. In addition, quite often it makes the code run faster. In particular,

# Logical expressions work not only on scalars, but also on vectors, matrices and, in general, on n-dimensional arrays.
# Subvectors/submatrices can be extracted using logical 0-1 arrays.

=== Irrelevant but useful example ===

typing <source enclose=none>a=1:5</source> in MATLAB command window creates a <math>1\times5</math> row-vector <source enclose=none>a</source> with values <math>[1\ 2\ 3\ 4\ 5]</math>. Logical expression <source enclose=none>ind=(a>3.5)</source> will create a so called logical vector <source enclose=none>ind</source> with values <math>[0\ 0\ 0\ 1\ 1]</math>, i.e. it is 1 if the according element is greater than 3.5 and 0 otherwise. Now, typing <source enclose=none>b=a(ind)</source> will generate a <math>2\times1</math> subvector <source enclose=none>b</source> with values <math>[4\ 5]</math>. You can also create some vectors or matrices with specific values changed: the command <source enclose=none>a(ind)=a(ind)*2</source> replace the last two values of the original vector <source enclose=none>a</source>. As a result, the vector <source enclose=none>a</source> becomes <math>[1 \ 2\ 3\ 8\ 10]</math>.

=== Slightly less irrelevant example ===

In some occasions you would like to modify the matrix of interest. Say, in some surveys “no answer” is coded as 999. Once you import the whole dataset in <source enclose=none>X</source>, you might want to replace these with, say, NaN. It can be done for the whole matrix of interest: <source enclose=none>X(X==999)=NaN</source>.

=== Relevant example ===

To demonstrate these capabilities in a more relevant environment, let’s run a very simple example. Assume that we have <math>T\times1</math> vector of returns <source enclose=none>r</source> and we want to

# Compute number of positive, negative and zero returns
# Compute means of positive and negative returns

The algorithm for this is quite straightforward:

# Find out the length of vector <source enclose=none>r</source>, T
# Initiate three counter variables, <source enclose=none>Tplus=0, Tzero=0, Tminus=0</source>, and vectors <source enclose=none>rplus=zeros(T,1), rminus=zeros(T,1)</source> (since we don’t know how many negative and positive returns we will observe
# Check whether r(i) is greater, smaller or equal to 0 for i=1
# If <source enclose=none>r(i)>0</source>, add 1 to Tplus, set <source enclose=none>rplus(Tplus)=r(i)</source>;
# Else if <source enclose=none>r(i)<0</source> add 1 to Tminus, set <source enclose=none>rminus(Tminus)=r(i)</source>;
# Else add 1 to Tzero
# Repeat 3-6 for <math>i=2,\ldots,T</math>
# Remove excessive zeros from <source enclose=none>rplus</source> and <source enclose=none>rminus</source>: <source enclose=none>rplus=rplus(1:Tplus);</source> 
<source enclose=none>rminus=rminus(1:Tminus);</source>
# Compute means of rminus and rplus. Number of positive, negative and zero returns are stored in <source enclose=none>Tplus,Tminus,Tzero</source>

MATLAB translation:

<source>T=size(r,1);
Tplus=0;Tminus=0;Tzero=0;
rplus=zeros(T,1);rminus=zeros(T,1);
for i=1:T
if r(i)>0
Tplus=Tplus+1;%increasing Tplus by one if return is positive
rplus(Tplus)=r(i);%storing positive return in the proper subvector
elseif r(i)<0
Tminus=Tminus+1;%increasing Tminus by one if return is negative
rminus(Tminus)=r(i);%storing negative return in the proper subvector
else
Tzero=Tzero+1;%increasing Tzero by one if return is neither positive nor negative
end
end
rplus=rplus(1:Tplus);%removing excessive zeros from a subvector of positive returns
rminus=rminus(1:Tminus);%removing excessive zeros from a subvector of negative returns
meanplus=mean(rplus);%computing mean of positive returns
meanminus=sum(rminus)/Tminus;%computing mean of negative returns</source>
Using MATLAB capabilities mentioned in this section, the algorithm can be reduced to:

# Construct a vector <source enclose=none>indplus</source> that has 1 for positive returns and 0 for negative returns
# Construct a vector <source enclose=none>indminus</source> that has 1 for negative returns and 0 for positive returns
# Assign to <source enclose=none>Tplus</source> a sum of elements of <source enclose=none>indplus</source>. This is the number of positive returns
# Assign to <source enclose=none>Tminus</source> a sum of elements of <source enclose=none>indminus</source>. This is the number of negative returns
# Compute <source enclose=none>Tzero</source> which is <source enclose=none>T-Tplus-Tminus</source>
# Construct a vector of positive returns <source enclose=none>rplus=r(indplus)</source> and compute its mean
# Construct a vector of negative returns <source enclose=none>rminus=r(indminus)</source> and compute its mean

MATLAB implementation:

<source> T=size(r,1);
indplus = r>0;%constructing an indicator vector with 1s if r(i)>0, 0 otherwise
indminus = r<0;%constructing an indicator vector with 1s if r(i)<0, 0 otherwise
Tplus=sum(indplus);%computing a number of positive returns
Tminus=sum(indminus);%computing a number of negative returns
Tzero=T-Tplus-Tminus;%computing a number of zero returns
rplus=r(indplus);%constructing a vector of positive returns
rminus=r(indminus);%constructing a vector of negative returns
meanplus=sum(rplus)/Tplus; %computing mean of positive returns
meanminus=mean(rminus); %computing mean of negative returns</source>
Or, a slightly shorter version of the same thing

<source> T=size(r,1);
rplus = r(r>0);%constructing a vector of positive returns
rminus = r(r<0);%%constructing a vector of negative returns
Tplus=size(rplus,1);%computing a number of positive returns
Tminus=size(indminus,l);%computing a number of negative returns
Tzero=T-Tplus-Tminus;%computing a number of zero returns
meanplus=sum(rplus)/Tplus; %computing mean of positive returns
meanminus=mean(rminus); %computing mean of negative returns</source>
A shorter code is less exposed to errors and easier to read (at least after some practice).

=Footnotes=

<references />

2013-10-11T12:52:20Z

AS: /* for ... end loop */

= Preliminaries =

Very often in your life you have to repeat the same operation many times (move your right and left leg sequentially while walking/running) or behave differently depending on external conditions (there is or there isn’t a bus on the bus stop). Quite often these two are combined together. Say, if there is a bus on the bus stop, then you run trying to catch it, otherwise walk or stop and enjoy the usual Manchester weather. The same is true for programming. Quite often you want to repeat the same operation many times, or you want to change the way you process your data depending on some conditions. We start with conditional statements. They execute different pieces of the code depending on whether <source enclose=none>condition</source> is true or false. There are several ways you can formulate it. The shortest

<source>if condition
statement1;
statement2;
...
end
</source>
, executes <source enclose=none>statement1;statement2,...</source> only if <source enclose=none>condition</source> is satisfied. <source enclose=none>condition</source> can be anything that generate non-zero or 0 (True or False), say <source enclose=none>i>0</source>, <source enclose=none>size(y,1)\sim=40</source>, or <source enclose=none>5-i</source>. The last condition is True always but for <source enclose=none>i=5</source>. A slightly longer version

<source> if condition
statement1;
statement2;
...
else
statement1a;
statement2a;
...
end
</source>
runs <source enclose=none>statement1;statement2;...</source>, if <source enclose=none>condition</source> is true and <source enclose=none>statement1a;statement2a;...</source> otherwise. The most general specification is

<source> if condition1
statement1;
statement2;
...
elseif condition2
statement1a;
statement2a;
...
...
...
elseif conditionN
statement1b;
statement2b;
...
else
statement1c;
statement2c;
...
end</source>
In this case, however, you have to ensure that <source enclose=none>condition1, condition2, …, conditionN</source> are mutually disjoint. As an example, you might think about different actions depending on your final grade. Condition1: <source enclose=none>grade<30</source>; Condition 2: <source enclose=none>(grade>=30)&&(grade<40)</source>; Condition 3: <source enclose=none>(grade>=40)&&(grade<50)</source>; etc.

MATLAB has two statements that create a loop. First, it is unconditional loop:

<source>for CounterVariable=[range of values]
statement1;
statement2;
...
end</source>
It repeats at most as many times as many elements it has in the <source enclose=none>[range of values]</source>. If the range of values is empty, this loop does not run. Say, if you define a range <source enclose=none>10:1</source>, MATLAB creates an empty range. Thus, this loop will not be executed. If you define a range <source enclose=none>1:3:10</source>, MATLAB creates a range of four values <source enclose=none>[1 4 7 10]</source>, and the loop runs four times. During the first iteration, <source enclose=none>CounterVariable=1</source>, during the second <source enclose=none>CounterVariable=4</source>, etc. After the end of the loop <source enclose=none>CounterVariable=10</source>. Please note, it is very unwise to modify the counter inside the loop. All modifications will disappear after the next iteration. Please also note, that the values in the range could be anything, including filenames or matrices from a cell vector. There are two commands that can modify the execution of the loop. <source enclose=none>continue</source> breaks the current ''iteration'' of the loop. Once it is executed, the loop continues skipping current iteration. <source enclose=none>break</source> stops the execution of the loop and your program continues after this point. These commands are used inside <source enclose=none>if</source> statements. For example, <source enclose=none>if CounterVariable == 10 continue;end</source> skips the loop iteration for <source enclose=none>CounterVariable = 10</source>.

Second, it is the conditional loop

<source>while condition
statement1;
statement2;
...
end</source>
This version of the loop executes statements as long as <source enclose=none>condition</source> is true. If <source enclose=none>condition</source> is always true, your loop runs forever.

== <source enclose=none>for ... end</source> loop ==

A standard application for <source enclose=none>for ... end</source> loop is the reconstruction of AR(p) series once AR(p) coefficients and the vector of error terms is known. <math>y_t=\phi_0+\sum_{i=1}^p \phi_i y_{t-i}+e_t.</math> For simplicity, we assume that <math>p=1</math>. Also, to be able to compute <math>y_1</math>, we need to provide <math>y_0</math>. Since we don’t know <math>y_0</math>,the best guess for <math>y_0</math> is <math>E(y_0)</math>. For stationary AR(1) process, that is for the case <math>|\phi_1|<1</math>, <math>E(y_0)=\phi_0/(1-\phi_1)</math>. Thus, knowing <math>y_0</math> and <math>e_t</math> for <math>t=1,\ldots,T</math>, we can reconstruct <math>y_t,\ t=1\ldots,T</math>:

<math>\begin{aligned}
y_1=&\phi_0+\phi_1 y_0+e_1\\
y_2=&\phi_0+\phi_1 y_1+e_2\\
&\ldots\\
y_t=&\phi_0+\phi_1 y_{t-1}+e_t\\
&\ldots\\
y_T=&\phi_0+\phi_1 y_{T-1}+e_T\end{aligned}</math>

Definitely, if you are patient enough and <math>T</math> is not very large, you can create your m file with <math>T</math> lines in it. However, once <math>T</math> is unknown, this approach would not work. Fortunately, there is a better alternative for this type of operations. All these computations can be summarized using the following algorithm, assuming <source enclose=none>e</source> is already created, say <source enclose=none>e=randn(1000,1)*2</source>:

# Find the length of a vector of error terms <source enclose=none>e</source>: <source enclose=none>T=size(e,1)</source>
# Initialize a vector <source enclose=none>y</source> of the same length as vector <source enclose=none>e</source>: <source enclose=none>y=zeros(T,1)</source>
# Compute <source enclose=none>y(1)=phi0+phi1*(phi0/(1-phi1))+e(1)</source>. Please remember, we assume that <math>y_0=E(y)=\phi_0/(1-\phi_1)</math>
# Compute <source enclose=none>y(i)=phi0+phi1*y(i-1)+e(i)</source> for <math>i=2</math>
# Repeat line 4 for <math>i=3,...,T</math>

When vector <source enclose=none>e</source> is known in advance, the MATLAB code is

<source> T=size(e,1);
y=zeros(T,1);
y0=phi0/(1-phi1);
y(1)=phi0+phi1*y0+e(1);
for i=2:T
y(i)=phi0+phi1*y(i-1)+e(i);
end</source>
However, if <source enclose=none>phi1=1</source>, <math>E(y_t)</math> is not constant, then the formula we use in the code does not work and will create either a series <source enclose=none>y</source> of <math>\pm\infty</math>, if <source enclose=none>phi0 \ne 0</source> or a series of not a numbers <source enclose=none>NaN</source>, if <source enclose=none>phi0=0</source><ref>There are two special numerical values in MATLAB. One is infinity <source enclose=none>Inf</source>, and another is not a number <source enclose=none>NaN</source>. A value of a variable becomes <source enclose=none>Inf</source> if the number is too big in absolute value (<math>\approx \pm 2e308</math>). Also, infinity is generated once you have expressions like <math>x/0</math>, where <math>x\ne0</math>. After that, infinity can only change a sign or become not a number. Not a number appears when there is an uncertainty of a kind of <math>0/0</math>, <math>\infty-\infty</math> and such. Any algebraic operations with <source enclose=none>NaN</source> result <source enclose=none>NaN</source>
</ref>.

== <source enclose=none>if else end</source> or <source enclose=none>if end</source> ==

To avoid these inconveniences, we have to consider separately two cases:

# AR(1) process is stationary, i.e. <math>|\phi_1|<1</math>
# AR(1) process is nonstationary, i.e. <math>|\phi_1|\ge 1</math>

For the latter, we have to acknowledge the fact that <math>E(y_t)=\mu_t</math>, i.e. unconditional expectation is a function of time. In this case we have to set <math>E(y_0)</math> to some value. A standard assumption for non-stationary series is to assume that <math>E(y_0)=0</math>.

The algorithm in this case would look like:

# Find a length of a vector of error terms <source enclose=none>e</source>: <source enclose=none>T=size(e,1)</source>
# Initialize a vector <source enclose=none>y</source> of the same length as vector <source enclose=none>e</source>: <source enclose=none>y=zeros(T,1)</source>
# Check whether <source enclose=none>abs(phi1)<1</source>. If this statement is true, then <source enclose=none>y0=phi0/(1-phi1)</source>. Else, <source enclose=none>y0=0</source>. Please remember, we set <math>y_0=E(y_0)</math>.
# Compute <source enclose=none>y(1)=phi0+phi1*y0+e(1)</source>.
# Compute <source enclose=none>y(i)=phi0+phi1*y(i-1)+e(i)</source> for <math>i=2</math>
# Repeat line 5 for <math>i=3,...,T</math>

Assuming vector <source enclose=none>e</source> is known in advance, the MATLAB code is

<source> T=size(e,1);
y=zeros(T,1);
if abs(phi1)<1
y0=phi0/(1-phi1);
else
y0=0;
end
y(1)=phi0+phi1*y0+e(1)
for i=2:T
y(i)=phi0+phi1*y(i-1)+e(i);
end</source>
If you don’t like the word else, you can skip it:

<source> T=size(e,1);
y=zeros(T,1);
y0=0;
if abs(phi1)<1
y0=phi0/(1-phi1);
end
y(1)=phi0+phi1*y0+e(1)
for i=2:T
y(i)=phi0+phi1*y(i-1)+e(i);
end</source>

== <source enclose=none>while end</source> loop ==

An alternative way of running the same code is to use a conditional loop (purely for demonstration purposes). Usually the conditional loop is used when the number of iterations is not known in advance.

# Find the length of a vector of error terms <source enclose=none>e</source>: <source enclose=none>T=size(e,1)</source>
# Initialize a vector <source enclose=none>y</source> of the same length as vector <source enclose=none>e</source>: <source enclose=none>y=zeros(T,1)</source>
# Check whether <source enclose=none>abs(phi1)<1</source>. If this statement is true, then <source enclose=none>y0=phi0/(1-phi1)</source>. Else, <source enclose=none>y0=0</source>. Please remember, we set <math>y_0=E(y_0)</math>.
# Compute <source enclose=none>y(1)=phi0+phi1*y0+e(1)</source>.
# Compute <source enclose=none>y(i)=phi0+phi1*y(i-1)+e(i)</source> for <math>i=2</math>
# Increase i by 1, i.e. <math>i=i+1</math> (please note, in programming this is an important statement),
# Repeat line 4 while <math>i<=T</math>

Assuming vector <source enclose=none>e</source> is known in advance, the MATLAB code is

<source> T=size(e,1);
y=zeros(T,1);
if abs(phi1)<1
y0=phi0/(1-phi1);
else
y0=0;
end
y(1)=phi0+phi1*y0+e(1)
i=2;
while i<=T
y(i)=phi0+phi1*y(i-1)+e(i);
i=i+1;
end</source>
== Imperfect substitutes of the above ==

MATLAB has two powerful tools that make programmer’s life much easier and utilization of loops/if less frequent. In addition, quite often it makes the code run faster. In particular,

# Logical expressions work not only on scalars, but also on vectors, matrices and, in general, on n-dimensional arrays.
# Subvectors/submatrices can be extracted using logical 0-1 arrays.

=== Irrelevant but useful example ===

typing <source enclose=none>a=1:5</source> in MATLAB command window creates a <math>1\times5</math> row-vector <source enclose=none>a</source> with values <math>[1\ 2\ 3\ 4\ 5]</math>. Logical expression <source enclose=none>ind=(a>3.5)</source> will create a so called logical vector <source enclose=none>ind</source> with values <math>[0\ 0\ 0\ 1\ 1]</math>, i.e. it is 1 if the according element is greater than 3.5 and 0 otherwise. Now, typing <source enclose=none>b=a(ind)</source> will generate a <math>2\times1</math> subvector <source enclose=none>b</source> with values <math>[4\ 5]</math>. You can also create some vectors or matrices with specific values changed: the command <source enclose=none>a(ind)=a(ind)*2</source> replace the last two values of the original vector <source enclose=none>a</source>. As a result, the vector <source enclose=none>a</source> becomes <math>[1 \ 2\ 3\ 8\ 10]</math>.

=== Slightly less irrelevant example ===

In some occasions you would like to modify the matrix of interest. Say, in some surveys “no answer” is coded as 999. Once you import the whole dataset in <source enclose=none>X</source>, you might want to replace these with, say, NaN. It can be done for the whole matrix of interest: <source enclose=none>X(X==999)=NaN</source>.

=== Relevant example ===

To demonstrate these capabilities in a more relevant environment, let’s run a very simple example. Assume that we have <math>T\times1</math> vector of returns <source enclose=none>r</source> and we want to

# Compute number of positive, negative and zero returns
# Compute means of positive and negative returns

The algorithm for this is quite straightforward:

# Find out the length of vector <source enclose=none>r</source>, T
# Initiate three counter variables, <source enclose=none>Tplus=0, Tzero=0, Tminus=0</source>, and vectors <source enclose=none>rplus=zeros(T,1), rminus=zeros(T,1)</source> (since we don’t know how many negative and positive returns we will observe
# Check whether r(i) is greater, smaller or equal to 0 for i=1
# If <source enclose=none>r(i)>0</source>, add 1 to Tplus, set <source enclose=none>rplus(Tplus)=r(i)</source>;
# Else if <source enclose=none>r(i)<0</source> add 1 to Tminus, set <source enclose=none>rminus(Tminus)=r(i)</source>;
# Else add 1 to Tzero
# Repeat 3-6 for <math>i=2,\ldots,T</math>
# Remove excessive zeros from <source enclose=none>rplus</source> and <source enclose=none>rminus</source>: <source enclose=none>rplus=rplus(1:Tplus);</source> 
<source enclose=none>rminus=rminus(1:Tminus);</source>
# Compute means of rminus and rplus. Number of positive, negative and zero returns are stored in <source enclose=none>Tplus,Tminus,Tzero</source>

MATLAB translation:

<source>T=size(r,1);
Tplus=0;Tminus=0;Tzero=0;
rplus=zeros(T,1);rminus=zeros(T,1);
for i=1:T
if r(i)>0
Tplus=Tplus+1;%increasing Tplus by one if return is positive
rplus(Tplus)=r(i);%storing positive return in the proper subvector
elseif r(i)<0
Tminus=Tminus+1;%increasing Tminus by one if return is negative
rminus(Tminus)=r(i);%storing negative return in the proper subvector
else
Tzero=Tzero+1;%increasing Tzero by one if return is neither positive nor negative
end
end
rplus=rplus(1:Tplus);%removing excessive zeros from a subvector of positive returns
rminus=rminus(1:Tminus);%removing excessive zeros from a subvector of negative returns
meanplus=mean(rplus);%computing mean of positive returns
meanminus=sum(rminus)/Tminus;%computing mean of negative returns</source>
Using MATLAB capabilities mentioned in this section, the algorithm can be reduced to:

# Construct a vector <source enclose=none>indplus</source> that has 1 for positive returns and 0 for negative returns
# Construct a vector <source enclose=none>indminus</source> that has 1 for negative returns and 0 for positive returns
# Assign to <source enclose=none>Tplus</source> a sum of elements of <source enclose=none>indplus</source>. This is the number of positive returns
# Assign to <source enclose=none>Tminus</source> a sum of elements of <source enclose=none>indminus</source>. This is the number of negative returns
# Compute <source enclose=none>Tzero</source> which is <source enclose=none>T-Tplus-Tminus</source>
# Construct a vector of positive returns <source enclose=none>rplus=r(indplus)</source> and compute its mean
# Construct a vector of negative returns <source enclose=none>rminus=r(indminus)</source> and compute its mean

MATLAB implementation:

<source> T=size(r,1);
indplus = r>0;%constructing an indicator vector with 1s if r(i)>0, 0 otherwise
indminus = r<0;%constructing an indicator vector with 1s if r(i)<0, 0 otherwise
Tplus=sum(indplus);%computing a number of positive returns
Tminus=sum(indminus);%computing a number of negative returns
Tzero=T-Tplus-Tminus;%computing a number of zero returns
rplus=r(indplus);%constructing a vector of positive returns
rminus=r(indminus);%constructing a vector of negative returns
meanplus=sum(rplus)/Tplus; %computing mean of positive returns
meanminus=mean(rminus); %computing mean of negative returns</source>
Or, a slightly shorter version of the same thing

<source> T=size(r,1);
rplus = r(r>0);%constructing a vector of positive returns
rminus = r(r<0);%%constructing a vector of negative returns
Tplus=size(rplus,1);%computing a number of positive returns
Tminus=size(indminus,l);%computing a number of negative returns
Tzero=T-Tplus-Tminus;%computing a number of zero returns
meanplus=sum(rplus)/Tplus; %computing mean of positive returns
meanminus=mean(rminus); %computing mean of negative returns</source>
A shorter code is less exposed to errors and easier to read (at least after some practice).

=Footnotes=

<references />