Andrea Bucci, Giulio Palomba and Marco Tedeschi 1.0 2026-05-18 Matrix Auto-Regressive Model C32 C40 Matrix AR(1) MAINWIN/Model/TSMulti pdfdoc:MAR.pdf examples Lists of Time Series Number of lags Model "with constant" "without constant" "without constant (demeaned time series)" Estimation method "Least Squares (LSE)" "Maximum Lkelihood (MLE)" Maximum number of 'flip-flop' iterations #Prepare data matrix SmplMean={} matrices DataMat=PrepareData(In,&SmplMean) if intercept==3 DataMat=DeMean(DataMat,SmplMean) endif #Set the bundle bdl=par(DataMat,lag,intercept) scalar m=bdl.m scalar n=bdl.n list allvars=empty loop i=1..n allvars+=In[i] endloop bdl.SampleMeans=SmplMean bdl.VarList=allvars bdl.maxiter=maxiter bdl.Vnames=varnames(allvars) #Simulate initial parameters: matrix Binit=InitBMat(&bdl) #Estimation: scalar nMat=1+2*lag matrices Theta=array(nMat) if est_method==1 string method="LSE" Estimate_LSE(Binit,&Theta,&bdl) elif est_method==2 string method="MLE" Estimate_MLE(Binit,&Theta,&bdl) endif matrix A=Theta[1] matrix B=Theta[2] matrix Mu=Theta[3] matrices Xhat=FittedTS(A,B,Mu,&bdl) matrices Ehat=Resids(&bdl) bdl.method=method matrix Σ=MakeSigma(&bdl) Loglik(&bdl) IC(&bdl) scalar rho=CheckStationarity(A,B) bdl.check=rho VCV(A,B,&bdl) PrintResults(A,B,Mu,&bdl) return 1 Row index Column index Forecast horizon (periods) Verbosity "NONE" "TABLE" "TABLE & PLOTS" string printout="PLOTS" if verb==1 printout="NONE" elif verb==2 printout="TABLE" endif scalar m = x.m scalar n = x.n strings Cnames=x.Vnames string warn="Generalized IRFs with orthogonal innovations cannot be computed." print " " if printout!="NONE" && printout!="TABLE" && printout!="PLOTS" printf "WARNING! ''%s'' is not a valid print command.\n",printout printf "Available print commands are ''NONE'', ''TABLE'' and ''PLOTS''. %s\n",warn print " " return 0 elif s_i>m && s_j<=n printf "WARNING! First argument greater than m=%g. %s\n",m,warn print " " return 0 elif s_j>n && s_i<=m printf "WARNING! Second argument greater than n=%g. %s\n",n,warn print " " return 0 elif s_i>m && s_j>n printf "WARNING! First and second arguments greater than m=%g and n=%g. %s.\n",m,n,warn print " " return 0 elif s_i<=0 || s_j<=0 printf "WARNING! First and second arguments greater must be positive. %s\n",warn print " " return 0 else matrix Sigma = x.Sigma matrix SigmaC = x.SigmaC matrix SigmaR = x.SigmaR matrix A = x.Theta[1] matrix B = x.Theta[2] scalar dim=m*n string method=x.method scalar H = horizon+1 matrix IRF = zeros(dim, H) strings Rnames=array(H) loop k=1..H Rnames[k]=sprintf("%g",k-1) endloop scalar pos=m*(s_j-1)+s_i loop k = 1..H if method == "LSE" IRF[,k] = (B^(k-1) ** A^(k-1))*Sigma[,pos] elif method == "MLE" IRF[,k] = (B^(k-1)*SigmaC[,s_j]) ** (A^(k-1)*SigmaR[,s_i]) endif endloop IRF=IRF' matrix out=IRF x.IRF=out string S=sprintf("%s---GIRFs(%g,%g)",Cnames[pos],s_i,s_j) if printout=="TABLE" || printout=="PLOTS" print "Shock-First Impulse Response Function with Orthogonal Innovations" printf "[%s]\n",S print " " rnameset(IRF,Rnames) cnameset(IRF,Cnames) print IRF strings names="t" names+=Cnames IRF=seq(1,H)'~IRF if printout=="PLOTS" cnameset(IRF, names) gpbuild IRFplot loop i=2..cols(IRF) gnuplot i 1 --matrix=IRF --with-lines --fit=none endloop end gpbuild gridplot IRFplot --fontsize=12 --height=800 --rows=n --title=@S --output=display endif endif return out endif Residual tests "All diagnostics" "Autocorrelation" "Heteroskedasticity" "Normality" Lags Multivariate tests string stringtest="ALL" if test==1 stringtest="AUTOCORR" elif test==2 stringtest="HETEROSK" elif test==3 stringtest="NORMALITY" endif scalar m=x.m scalar n=x.n scalar T=x.T scalar p=x.p+1 matrices Ehat=x.Resids matrix E=zeros(p-1,m*n) loop t=p..T E|=vec(Ehat[t])' endloop list R=E matrix out={} var 0 R --silent if multivariate print "Multivariate Diagnostic tests:" if stringtest=="AUTOCORR" modtest lag --autocorr out=seq(1,lag)'~$test~$pvalue print " " elif stringtest=="HETEROSK" modtest lag --arch out=seq(1,lag)'~$test~$pvalue print " " elif stringtest=="NORMAL" modtest --normality out=NA~$test~$pvalue print " " elif stringtest=="ALL" modtest lag --autocorr out|=seq(1,lag)'~$test~$pvalue print " " modtest lag --arch out|=seq(1,lag)'~$test~$pvalue print " " modtest --normality --quiet out|=NA~$test~$pvalue endif else print "Univariate Diagnostic tests on vec(E):" if stringtest=="AUTOCORR" modtest lag --autocorr --univariate out|=$test~$pvalue print " " elif stringtest=="HETEROSK" modtest lag --arch --quiet --univariate out|=$test~$pvalue print " " elif stringtest=="NORMAL" print "Distributional shape test:" print " " loop foreach i R normtest R[i] --jbera out|=$test~$pvalue endloop print " " elif stringtest=="ALL" modtest lag --autocorr --univariate out|=$test~$pvalue print " " modtest lag --arch --quiet --univariate out|=$test~$pvalue print " " print "Distributional shape test:" print " " loop foreach i R normtest R[i] --jbera out|=$test~$pvalue endloop print " " endif endif x.Diagnostics=out return out # This function sets the bundle # In: (1) a matrix with data in time series setting: t-th row (T rows) --> put all variables for each country (m variables for n countries) # (2) MAR order (n. of lags) # (3) int: if 1 the model has a constant, if 2 the model has not a constant, if 3 if 2 the model has not a constant and time series are demeaned # Out: the bundle bundle ret ret["AllDataMat"]=DataSet ret["m"]=rows(DataSet[1]) ret["n"]=cols(DataSet[1]) ret["T"]=$nobs ret["withconst"]=intercept ret["iterations"]=1 ret["p"]=lag string s1=sprintf("%d",$obsdate[1+lag]) string s2=sprintf("%d",$obsdate[ret["T"]]) string ret["startdate"]=sprintf("%s-%s-%s",s1[1:4],s1[5:6],s1[7:8]) string ret["enddate"]=sprintf("%s-%s-%s",s2[1:4],s2[5:6],s2[7:8]) return ret # This function transforms the lists of variables in an array of m x n matrices scalar n=nelem(In) scalar m=nelem(In[1]) scalar T=$nobs matrices Data=empty matrix SampleMean=0 matrix tmp={} loop t=1..T loop i=1..n loop j=1..m tmp~=In[i][j][t] endloop endloop SampleMean+=tmp Data+=mshape(tmp,m,n) tmp={} endloop SM=mshape(SampleMean,m,n)/T return Data scalar T=nelem(M) matrices out=empty loop i=1..T out+=M[i]-Mean endloop return out # This function standardises a matrix by using its Frobenius norm # In: a matrix # Out: standardised matrix with the same dimensions as In return In/sqrt(sumr(sumc(In.^2))) # This function performs the LSE estimation via the the so-called `flick-flack' methodology # In: (1) the initial matrix B # (2) an address with the parameters vector # (3) the bundle scalar m=x.m scalar n=x.n scalar T=x.T scalar p=x.p+1 scalar K=x.maxiter matrices X=x.AllDataMat scalar i=x.iterations scalar threshold=10^(-9) matrix Mu=0 if x.withconst==1 matrix Mu=x.SampleMeans endif matrix B=StartB matrix oldMu=Mu matrix oldA=zeros(m,m) matrix oldB=zeros(n,n) loop while i<K matrix tmpA1=zeros(m,m) matrix tmpA2=zeros(m,m) matrix tmpB1=zeros(n,n) matrix tmpB2=zeros(n,n) loop t=p..T tmpA1+=(X[t]-Mu)*B*X[t-1]' tmpA2+=X[t-1]*B'B*X[t-1]' endloop A=ScaleMatrix(tmpA1/tmpA2) loop t=p..T tmpB1+=(X[t]-Mu)'*A*X[t-1] tmpB2+=X[t-1]'*A'A*X[t-1] endloop B=tmpB1/tmpB2 if x.withconst==1 loop t=p..T Mu+=X[t]-A*X[t-1]*B endloop Mu=Mu/(T-p+1) endif scalar mm=abs(sumr(sumc(Mu-oldMu))) scalar aa=abs(sumr(sumc(A-oldA))) scalar bb=abs(sumr(sumc(B-oldB))) if aa<threshold && bb<threshold x.iterations=i i=K else oldMu=Mu oldA=A oldB=B i++ x.iterations=K endif endloop Params[1]=A Params[2]=B Params[3]=Mu matrices x.Theta=Params #compute Σr and Σc scalar divR=n*(T-x.p) scalar divC=m*(T-x.p) FittedTS(A,B,Mu,&x) matrices Ehat=Resids(&x) matrix SigmaC = InitSigmaC(&x) matrix SigmaR = empty loop i=1..10 matrix tmpSigmaR=zeros(m,m) matrix tmpSigmaC=zeros(n,n) loop t=p..T tmpSigmaR+=Ehat[t]/SigmaC*Ehat[t]' endloop SigmaR=ScaleMatrix(tmpSigmaR/divR) loop t=p..T tmpSigmaC+=Ehat[t]'/SigmaR*Ehat[t] endloop SigmaC=tmpSigmaC/divC endloop x.SigmaR=SigmaR x.SigmaC=SigmaC x.SigmaKron=SigmaC**SigmaR # This function performs the MLE estimation via the the so-called `flick-flack' methodology (sec 3.3 paper) # In: (1) the initial matrix B # (2) an address with the parameters vector # (3) the bundle #note: this is slightly different from Iterated computations given different Σ scalar K=x.maxiter scalar m = x.m scalar n = x.n scalar T = x.T scalar p=x.p+1 matrices X=x.AllDataMat matrix Mu=0 if x.withconst==1 matrix Mu=x.SampleMeans endif matrix A=I(m)/sqrt(m) matrix B=StartB scalar divR=n*(T-x.p) scalar divC=m*(T-x.p) FittedTS(A,B,Mu,&x) Resids(&x) matrix SigmaC = InitSigmaC(&x) matrix SigmaR = empty matrices Ehat=empty scalar i=x.iterations matrix oldMu=Mu matrix oldA=zeros(m,m) matrix oldB=zeros(n,n) scalar threshold=10^(-9) loop while i<K matrix tmpA1=zeros(m,m) matrix tmpA2=zeros(m,m) matrix tmpB1=zeros(n,n) matrix tmpB2=zeros(n,n) matrix tmpSigmaR=zeros(m,m) matrix tmpSigmaC=zeros(n,n) loop t=p..T tmpA1+=(X[t]-Mu)/SigmaC*B*X[t-1]' tmpA2+=X[t-1]*qform(B',inv(SigmaC))*X[t-1]' endloop A=ScaleMatrix(tmpA1/tmpA2) FittedTS(A,B,Mu,&x) Ehat = Resids(&x) loop t=p..T tmpSigmaR+=Ehat[t]/SigmaC*Ehat[t]' endloop SigmaR=ScaleMatrix(tmpSigmaR/divR) loop t=p..T tmpB1+=(X[t]-Mu)'/SigmaR*A*X[t-1] tmpB2+=X[t-1]'*qform(A',inv(SigmaR))*X[t-1] endloop B=tmpB1/tmpB2 if x.withconst==1 loop t=p..T Mu+=X[t]-A*X[t-1]*B endloop Mu=Mu/(T-p+1) endif FittedTS(A,B,Mu,&x) Ehat = Resids(&x) loop t=p..T tmpSigmaC+=Ehat[t]'/SigmaR*Ehat[t] endloop SigmaC=tmpSigmaC/divC scalar mm=abs(sumr(sumc(Mu-oldMu))) scalar aa=abs(sumr(sumc(A-oldA))) scalar bb=abs(sumr(sumc(B-oldB))) if aa<threshold && bb<threshold x.iterations=i i=K else oldMu=Mu oldA=A oldB=B i++ x.iterations=K endif endloop x.SigmaR=SigmaR x.SigmaC=SigmaC x.SigmaKron=SigmaC**SigmaR Params[1]=A Params[2]=B Params[3]=Mu matrices x.Theta=Params # This function computes the loglikelihhod # In: the bundle # Out: the loglikelihhod matrices Ehat=x.Resids matrix SigmaR=x.SigmaR matrix SigmaC=x.SigmaC scalar m=x.m scalar n=x.n scalar T=x.T scalar p=x.p+1 scalar divR=n*(T-x.p) scalar divC=m*(T-x.p) scalar LLK = 0 loop t =p..T LLK += tr((SigmaR\Ehat[t])*(SigmaC\Ehat[t]')) endloop LLK =-(LLK+divC*log(det(SigmaC))+divR*log(det(SigmaR))) x.LLK = LLK return LLK # This function computes the information criteria # In: the bundle # Out: vector whose components are AIC, BIC and HQC scalar m=x.m scalar n=x.n scalar T=x.T-x.p scalar k=m^2+n^2 if x.withconst==1 k=k+m*n endif scalar tmp=-2*x.LLK scalar AIC=tmp+2*k scalar BIC=tmp+k*log(T) scalar HQC=tmp+2*k*log(log(T)) x.AIC=AIC x.BIC=BIC x.HQC=HQC matrix out=AIC|BIC|HQC return out # This function simulates data given # In: (1) an initial Data Matrix (X0) # (2) a given parameter matrix A # (3) a given parameter matrix B # (4) a scalar indicating the sample size # (5) a scalar indicating the lags # Out: list of simulated time series scalar pp = p+1 scalar m = cols(A) scalar n = cols(B) scalar k = m*n matrices Y=empty matrix y = zeros(T, m*n) Y+=X0 matrix y[1,] = vec(X0)' loop t=pp..T Y+=A*Y[t-1]*B'+mnormal(m,n) y[t,]= vec(Y[t])' endloop strings V = array(k) loop i=1..k V[i] = sprintf("V%d", i) endloop cnameset(y,V) list L = mat2list(y) return L # This function computes an array of matrices with fitted data # In: (1) estimated parameter matrix A # (2) estimated parameter matrix B # (3) estimated parameter matrix Mu # (3) the bundle # Out: T matrices of dimension m x n containing fitted X^=Mu+A^*X(t-1)*B^' scalar p=x.p+1 matrices out=zeros(x.m,x.n) scalar T=x.T matrices Data=x.AllDataMat loop t=p..T out+=parMu+parA*Data[t-1]*parB' endloop x.Fitted=out return out # This function computes an array of matrices with residuals # In: the bundle # Out: T matrices of dimension m x n containing E^=Y(t)-A^*Y(t-1)*B^' matrices out=zeros(x.m,x.n) scalar T=x.T scalar p=x.p+1 matrices Data=x.AllDataMat matrices Fit=x.Fitted loop t=p..T out+=Data[t]-Fit[t] endloop x.Resids=out return out # This function sets the initial value of the column covariance matrix of residuals Σc # In: the bundle # Out: initial matrix Σc with dimension n x n scalar T=x.T scalar m=x.m scalar n=x.n scalar p=x.p+1 matrix out=zeros(n,n) matrices Ehat=x.Resids loop t=p..T out+=Ehat[t]'Ehat[t] endloop out=out/((T-x.p)*m) x.InitSigma=out return out # This function returns the parameter covariance matrix # In: (1) estimated parameter matrix A # (2) estimated parameter matrix B # (3) estimated parameter matrix Mu # (4) the bundle # Out: m^2+n^2 matrix scalar m=x.m scalar n=x.n scalar T=x.T scalar p=x.p+1 matrices Data=x.AllDataMat scalar dim=m^2+n^2 matrix γ=vec(parA)|zeros(n^2,1) if x.withconst==1 dim=dim+m*n γ=zeros(m*n,1)|vec(parA)|zeros(n^2,1) endif matrix V=zeros(dim,dim) matrix V1=zeros(dim,dim) matrix W={} matrix W1={} string method=x.method if method=="LSE" matrix Sigma=x.Sigma loop t=p..T W=(Data[t-1]*parB')**I(m)|I(n)**(Data[t-1]'parA') if x.withconst==1 W=I(m*n)|W endif V1+=qform(W,Sigma) V+=W*W' endloop matrix H=V+(T-x.p)*γ*γ' matrix Xi=qform(inv(H),V1) elif method=="MLE" matrix Sigma=x.SigmaKron loop t=p..T W=(Data[t-1]*parB')**I(m)|I(n)**(Data[t-1]'parA') if x.withconst==1 W=I(m*n)|W endif V+=qform(W,inv(Sigma)) endloop #matrix H=V+(T-x.p)*γ*γ' #matrix Xi=qform(inv(H),V) #this is Xi3 in the original paper V=V/(T-x.p) matrix H=(V+γ*γ') matrix Xi=qform(inv(H),V) Xi=Xi/(T-x.p) endif # Reshaping mechanism for the matrix B parameter standard errors scalar until=m^2 if x.withconst==1 until=until+m*n endif matrix id=seq(1,until)'|(until+ReshapeMatrix(n)) Xi=Xi[id,id] matrix x.VCV=Xi matrix x.SE=sqrt(diag(Xi)) return Xi # This function prints the estimation output # In: (1) estimated parameter matrix A # (2) estimated parameter matrix B # (2) estimated parameter matrix Mu # (4) the bundle # Out: matrix of dimension [mxn+p(m^2+n^2)] x 4 with the estimation output (without the intercept the dimension is p(m^2+n^2) x 4) scalar m=x.m scalar n=x.n scalar T=x.T-x.p scalar df=x.T-x.p*(m+n) scalar max_iter=x.maxiter scalar n_iter=x.iterations scalar dim=m^2+n^2 scalar c_dim=0 strings VarNames=empty string ini=x.startdate string fin=x.enddate string method=x["method"] print " " printf "MAR(1), using observations %s:%s (T = %g).\n",ini,fin,T if x.withconst==3 printf "Each data matrix at time t has dimension %gx%g (demeaned time series).\n",m,n else printf "Each data matrix at time t has dimension %gx%g.\n",m,n endif printf "Estimation method: %s.\n",method if n_iter<max_iter printf "Convergence achieved after %g iterations.\n",n_iter else printf "WARNING! Maximum number of %g iterations reached. Convergence not achieved.\n",max_iter endif print " " if x.withconst==1 c_dim=m*n matrix pars=vec(parMu)|vec(parA)|vec(parB) matrix se=x.SE matrix EstimationOutput=pars~se matrix outMu=EstimationOutput[1:c_dim,] loop i=1..n loop j=1..m VarNames+=sprintf("const[%g,%g]",j,i) endloop endloop print "Intercepts:" modprint outMu VarNames else matrix pars=vec(parA)|vec(parB) matrix se=x.SE matrix EstimationOutput=pars~se endif VarNames=empty loop i=1..m loop j=1..m VarNames+=sprintf("A[%g,%g]",j,i) endloop endloop matrix outA=EstimationOutput[c_dim+1:c_dim+m^2,] print "Matrix A:" modprint outA VarNames VarNames=empty loop i=1..n loop j=1..n VarNames+=sprintf("B[%g,%g]",j,i) endloop endloop matrix outB=EstimationOutput[c_dim+m^2+1:,] scalar rho=x.check print "Matrix B:" modprint outB VarNames if x.LLK>0 printf "Log-likelihood%14.7g Akaike criterion%10.7g\n",x.LLK,x.AIC else printf "Log-likelihood%13.7g Akaike criterion%10.7g\n",x.LLK,x.AIC endif if x.BIC>0 printf "Schwarz criterion%10.7g Hannan-Quinn%14.7g\n\n",x.BIC,x.HQC else printf "Schwarz criterion%9.7g Hannan-Quinn%14.7g\n\n",x.BIC,x.HQC endif print "Check for time series stationarity:" printf " ρ(A)ρ(B)=%.5f\n",rho if rho>=1 printf "WARNING! The index is more than one, the estimated MAR(%g) model is not stationary and causal.\n",x.p else printf "The estimated MAR(%g) model is stationary and causal.\n",x.p endif print " " if x.method=="LSE" matrix Sigma=x.Sigma scalar lambda=eigensym(Sigma)[1] if lambda<=0 print "WARNING! The estimated residual covariance matrix is not positive (semi)definite." endif elif x.method=="MLE" matrix SigmaK=x.SigmaKron matrix Sigmar=x.SigmaR matrix Sigmac=x.SigmaC scalar lambda=eigensym(SigmaK)[1] scalar lambdar=eigensym(Sigmar)[1] scalar lambdac=eigensym(Sigmac)[1] if lambda<=0 print "WARNING! The estimated residual covariance Σ=Σc**Σr is not positive (semi)definite." endif if lambdar<=0 print "WARNING! The estimated residual covariance Σr is not positive (semi)definite." endif if lambdac<=0 print "WARNING! The estimated residual covariance Σc is not positive (semi)definite." endif endif print " " return EstimationOutput # This function computes the causality condition (see Proposition 1 in the paper). # Specifically, it return the product of the spectral radii of estimated matrices A and B matrix out1=abs(eigen(parA)) matrix out2=abs(eigen(parB)) scalar out=max(out1)*max(out2) return out # This function computes the residual covariance matrix Σ=vec(E)*vec(E)/(T-p) scalar dim=x.m*x.n scalar T=x.T scalar p=x.p+1 matrices Ehat=x.Resids matrix e={} matrix out=zeros(dim,dim) loop t=p..T e=vec(Ehat[t]) out+=e*e' endloop x.Sigma=out/(T-x.p) return out/(T-x.p) # Reshaping mechanism for the matrix B parameter stanard errors scalar n=dim matrix out=vec((ones(n,1)*seq(0,n^2-1,n))'+ones(n,1)*seq(1,n)) return out # This function sets an initialisation for the parameter matrix B scalar m=x.m scalar n=x.n matrix Start=mrandgen(u,0.05,0.1,n,n) # Sometimes the rank of the square matrix B is not full. This is a correction mechanism for this issue. if rank(Start)==n Start+=mrandgen(u,0.01,0.05,n,n) endif return Start # This function reshapes a square matrix # In: square matrix B**A' (dim: mn x mn) # Out: matrix vec(A)vec(B)' (dim. m^2 x n^2) [AT THE MOMENT] matrix tmp=zeros(m^2,n^2) scalar col=1 matrix pr=empty matrix pc=empty loop i=1..n pr=seq((i-1)*m+1,i*m) loop j=1..n pc=seq((j-1)*m+1,j*m) tmp[,col]=mshape(X[pr,pc],m^2,1) col++ endloop endloop return tmp # This function carries out the specification test (see section 4.2 in the paper) # In: if printout=1 prints the test # note1: theoretically, it is valid only for the 'PROJ' estimator (not for LSE and MLE!!!) # note2: to check... scalar m=x.m scalar n=x.n scalar T=x.T scalar p=x.p scalar k=m*n matrix alpha=vec(x.Theta[1]) matrix beta=vec(x.Theta[2]) list L=x.VarList var 1 L --nc --silent matrix Phi_hat=$coeff' matrix Σ=$sigma matrix y={L} matrix Gamma0=mcov(y,T-p) matrix id=mshape(seq(1,k^2),k,k)' id=vec(Kron2vec(id,m,n)) matrix Ξ=(inv(Gamma0)**Σ)[id,id] matrix Phi_tilde=Kron2vec(Phi_hat,m,n) matrix vAvB=alpha*beta' beta=vec(ScaleMatrix(x.Theta[2])) matrix Dhat=Phi_tilde-vAvB matrix P=(I(rows(beta))-beta*beta')**(I(rows(alpha))-alpha*alpha') matrix V=ginv(P*Ξ*P) scalar test=T*vec(Dhat)'V*vec(Dhat) scalar df=(m^2-1)*(n^2-1) scalar pval=pvalue(X,df,test) matrix out=test~df~pval matrix x.SpecifTest=out if printout print " " print "Specification test" print "Null hypothesis: MAR(1) is preferable to VAR(1) on vecs" printf "Test statistic: %f\n",test printf "with p-value = P(Chi-Square(%g) > %f) = %g\n",df,test,pval print " " endif return out set verbose off include MAR.gfn #Opening the file containing data open Datavec.gdt --frompkg=MAR #Setting a seed for random number generation set seed 15061973 #Choosing the maximum number of 'flip-flop' iterations scalar maxiter = 100 #Number of lags scalar p=1 #constant matrix in the MAR(1) scalar intercept=2 #Defining the bundle bundle bdl # Preparation of n lists (in columns) containing m variables (in rows) list S1=V1 V2 list S2=V3 V4 list S3=V5 V6 list S4=V7 V8 #Creating the array of lists required by the public function lists X=defarray(S1,S2,S3,S4) # LSE estimation MAR(X,p,intercept,1,maxiter,&bdl) # MLE estimation MAR(X,p,intercept,2,maxiter,&bdl) # Multivariate diagnostics MAR_Diagnostics(,4,1,&bdl) # Multivariate diagnostics MAR_Diagnostics(,4,0,&bdl) # GIRFs MAR_IRFs(1,2,10,,&bdl) print "-----------------------------------------------------------------" # Print the contents of the bundle print bdl #---------------------------------------------------------------------------