' Simulation parameter

!N=1000		'Number of draws
!T=1000		'Sample size
!wind=100		'Window 
!arl=12  			'Max AR comparison
!tn=!T+!wind

rndseed 123   	 'Seed

create u	!tn				 'Creating workspace

vector(!N) icbic= 0
vector(!N) icaic= 0
vector(!N) ichq= 0
matrix(!N,!arl) dicbic=0
matrix(!N,!arl) dicaic= 0
matrix(!N,!arl) dichq= 0

smpl 1 !tn   				 'Range to work
for !i = 1 to !N

vector(!arl) selbic= 0
vector(!arl) selaic= 0
vector(!arl) selhq= 0

'DGP
'**************************************************
' y is AR(3)	
		'Parameters:
		!c=2
		!phi1=0.3
		!phi2=0.3
		!phi3=0.1

genr y= !c + nrnd

smpl 4 !tn
genr y= !c + !phi1*y(-1) + !phi2*y(-2) + !phi3*y(-3) + nrnd

smpl !wind+1 !tn

for !j = 1 to !arl 			
equation eq{!j}.ls y c y(-1 to -!j)
selbic(!j)=eq{!j}.@schwarz 
selaic(!j)=eq{!j}.@aic
selhq(!j)=eq{!j}.@hq

next

vector posbic=@ranks(selbic)
vector posaic=@ranks(selaic)
vector poshq=@ranks(selhq)

for !l = 1 to !arl
if posbic(!l)=1 then
icbic(!i)=!l
endif
if posaic(!l)=1 then
icaic(!i)=!l
endif
if poshq(!l)=1 then
ichq(!i)=!l
endif

if posbic(!l)=1 then
dicbic(!i,!l)=1
endif
if posaic(!l)=1 then
dicaic(!i,!l)=1
endif
if poshq(!l)=1 then
dichq(!i,!l)=1
endif

next

'**************************************************
next

range 1 !N  'Simulation results range
smpl 1 !N

mtos(icbic,bic)
mtos(icaic,aic)
mtos(ichq,hq)

vector bicfreq = @cmeans(dicbic)
vector aicfreq = @cmeans(dicaic)
vector hqfreq = @cmeans(dichq)

graph hist_ics.distplot(scale="relfreq") bic aic hq
hist_ics.axis(b) range(1,10)
hist_ics.axis(l) range(0,1000)

graph dot_ics.dot bic aic hq

scalar tsize=!arl+2
table(tsize,4) Results
Results(1,1)="Lag length"
Results(1,2)="BIC"
Results(1,3)="AIC"
Results(1,4)="HQ"

for !i=3 to tsize
Results(!i,1)=!i-2
Results(!i,2)=bicfreq(!i-2)
Results(!i,3)=aicfreq(!i-2)
Results(!i,4)=hqfreq(!i-2)
next
Results.setlines(2) d
Results.setformat(@all) ft.2  

delete dicaic dicbic dichq icaic icbic ichq posaic posbic poshq selaic selbic selhq tsize 

show hist_ics
show dot_ics
Results.display