Analyza casoveho radu poctu novozaregistrovanych aut v USA, mesacne data

SetDirectory["F:\\Analyza.Casovych.Radov"]

SetOptions[ListPlot,PlotRange->All,AxesLabel->{"time [Year]","pocet * 1000"},BaseStyle->{11,FontFamily->"Arial"},Joined->True];

nacitame mesacne udaje o pocte zaregistrovnych novych aut v USA ( v tisicoch) za roky 1982 az 1988.

data = ReadList["auta_mes.dat",Number];

n=Length[data]

Save["auta.def",data]

ZAKLADNA ANALYZA CASOVEHO RADU

vypocet vyberovej strednej hodnoty, smerodajnej odchylky a rozptylu casoveho radu

μ=Mean[data]//N
σ=StandardDeviation[data]//N
σ2=Variance[data]//N

popis casovej osi pri zobrazeni casoveho radu mesacnych prirastkov novoregistrovanych aut

popis=Table[{13+i*12,1983+i},{i,0,5}];

gd=ListPlot[data,Ticks->{popis,Automatic},AxesOrigin->{0,0}]

vypocet casoveho radu rocnych prirastkov novoreg. aut.
popis casovej osi (s krokom 2 roky) pri zobrazeni casoveho radu.
zobrazenie dat. ulozenie grafu do premennej.

dat=Partition[data,12];

rp - pocty predanych aut za rok v tisicoch

rp=Table[(∑_(j=1)^12)dat[[i,j]],{i, n/12}]

popisr=Table[{1+i*2,1983+i*2},{i,0,5}];

gdr=ListPlot[rp,Ticks->{popisr,Automatic},PlotStyle->Hue[0.05],AxesOrigin->{0,0}]

vycislime priemerne prirastky aut v jednotlivych mesiacoch (za cele sledovane obdobie)

pocetr=Length[dat]

mesiac=Table[i,{i,1,12}];

pmp=N[Table[(1/pocetr)(∑_(i=1)^pocetr)dat[[i,j]],{j,12}]]

matica=Transpose[{mesiac,pmp}];

TableForm[matica,TableHeadings->{None,{"mesiac","priemer"}}]

ListPlot[pmp,AxesLabel->{"month","pocet * 1000"},AxesOrigin->{0,700}]

z grafu je zrejme, ze najmenej aut sa predava vo februari

URCENIE SYSTEMATICKYCH ZLOZIEK CASOVEHO RADU

<<"LinearRegression`"

Needs["NonlinearRegression`"]

x=data;

Trend

a) vycislenie a zobrazenie prvych diferencii

V pripade, ze su prve diferencie priblizne nulove, pouzijeme konst. trend, ak su prve diferencie priblizne konst. (ale nenulove), pouzijeme linearny trend.

tv1=Drop[x,1]-Drop[x,-1];

Variance[tv1]//N
ListPlot[tv1,AxesLabel->{"time [month]","flow [m^3/s]"},Joined->False]

b) vycislime a zobrazime druhe diferenicie casoveho radu

tv2=Drop[x,2]-2*Drop[Drop[x,-1],1]+Drop[x,-2];
Variance[tv2]//N
ListPlot[tv2,AxesLabel->{"time [month]","flow [m^3/s]"},Joined->False]

c) vycislime a zobrazime podiely susednych clenov radu

v pirpade, ze su podiely

tv3=Table[If[(x[[t]]≠0)∧(x[[t+1]]≠0),(1/x[[t+1]])/(1/x[[t]]),0],{t,n-1}];
Variance[tv3]//N
ListPlot[tv3,AxesLabel->{"time [month]","flow [m^3/s]"},Joined->False]

vacsina podielov lezi v intervale (0;1.5), preto mozno pouzit exponencialny trend

d) vycislime a zobrazime podiely na testovanie vhodnosti logistickeho trendu

tv4=Table[If[(x[[t]]≠0)∧(x[[t+1]]≠0)∧(x[[t+2]]≠0)∧(1/x[[t+1]]-1/x[[t]]≠0),(1/x[[t+2]]-1/x[[t+1]])/(1/x[[t+1]]-1/x[[t]]),0],{t,n-2}];
Variance[tv4]//N
ListPlot[tv4,AxesLabel->{"time [month]","flow [m^3/s]"},Joined->False]

2. Vypocet trendu

Exponencialny trend

rege=NonlinearRegress[x,a*b^t,{a,b},t,RegressionReport->{BestFit,ParameterCITable,EstimatedVariance,ANOVATable}]

etrend=BestFit/.rege

gte=Plot[etrend,{t,0,n},PlotStyle->{Hue[0.6],Thickness[0.005]}];
Show[gd,gte]

rege=NonlinearRegress[x,a*b^t,{a,b},t,RegressionReport->{FitResiduals}];
rezte=FitResiduals/.rege;

μ=Mean[rezte]//Chop
σ=StandardDeviation[rezte]
σ2=Variance[rezte]

ListPlot[rezte,PlotStyle->{Hue[0.6],Thickness[0.005]}]

Logisticky trend

regl=NonlinearRegress[x,γ/(1+a*b^t),{{a,-1},{b,0.2},{γ,50}},t,RegressionReport->{BestFit,ParameterCITable,EstimatedVariance,ANOVATable}]

ltrend=BestFit/.regl

gtl=Plot[ltrend,{t,0,n},PlotStyle->{Hue[0.6],Thickness[0.005]}];
Show[gd,gtl]

regl=NonlinearRegress[x,γ/(1+a*b^t),{{a,-1},{b,0.2},{γ,50}},t,RegressionReport->{FitResiduals}];

reztl=FitResiduals/.regl;

vypocet vyberovej str. hodn., smerod. odch. a roptylu casoveho radu reziui po logistickom trende

Save["reztl.dat",reztl]

15. april  2008

reztl=ReadList["reztlauta.dat"][[1]];

<<"LinearRegression`"

Needs["NonlinearRegression`"]

<<"timeseri\\datasmoo.m"
<<"timeseri\\timeseri.m"
<<"timeseri\\userfunc.m"

Analyzujeme : reztl

Odhad sezonnej zlozky pomocou korelogramu :

kf = CorrelationFunction[reztl, 60];
myplotcorr1[kf,2/Sqrt[n],AxesLabel->{k,(Subscript[ρ, k])},Ticks->{{12,24,36,48,60},Automatic},PlotLabel-> "AutoCorrelation Function",PlotRange->All]

Korelogram ma pravidelne extremy. Teda rezidua obsahuju sezonnu zlozku.

Urcenie sezonnej zlozky pomocou matematickych funkcii :

Regress[reztl,{Sin[2*Pi*t/12],Cos[2*Pi*t/12]},t,IncludeConstant->False]

regmf=Regress[reztl,{Sin[2*Pi*t/12],Cos[2*Pi*t/12]},t,IncludeConstant->False,RegressionReport->{BestFit,FitResiduals}];

sz=BestFit/.regmf

rezszmf=FitResiduals/.regmf

ListPlot[rezszmf,Ticks->{popis,Automatic},Joined->True,AspectRatio->1/3]

μ=Mean[rezszmf]//Chop
σ=StandardDeviation[rezszmf]
σ2=Variance[rezszmf]

Cyklicka zlozka

spec=Spectrum[rezszmf];
gs=plotspectrum[spec,AxesLabel->{"ω","f(ω)"},Joined->True]
dv=Transpose[{spec,Table[i,{i,0,n-1}]}];
Take[Sort[dv],-20]

Vyznamne frekvencie urcime pomocou Fischerovho perioicity :

M=Round[n/2];
α=0.05;
mF=M;
Y[j_]:=dv[[j,1]]/Sum[dv[[i,1]],{i,M}]
sY1=Take[Sort[Table[{Y[j],j-1},{j,M}]],{-1}]
1-(α/mF)^(1/(mF-1))

dlzkaPeriody = n/28.

(Subscript[ω, 28])=dlzkaPeriody;

vf=(2*Pi *28/n);
funkcie=Flatten[Table[{Sin[vf[[i]]*t],Cos[vf[[i]]*t]},{i,Length[vf]}]];
Regress[rezszmf,funkcie,t,IncludeConstant->False]

regcz=Regress[reztl,funkcie,t,IncludeConstant->False,RegressionReport->{BestFit,FitResiduals}];

cz=BestFit/.regcz

rez=FitResiduals/.regcz

ListPlot[rez,Joined->True,kf = CorrelationFunction[reztl, 60];
myplotcorr1[kf,2/Sqrt[n],AxesLabel->{k,(Subscript[ρ, k])},Ticks->{{12,24,36,48,60},Automatic},PlotLabel-> "AutoCorrelation Function",PlotRange->All]

Syst

corrx= CorrelationFunction[rez, 60];
myplotcorr1[corrx,2/Sqrt[n],AxesLabel->{k,(Subscript[ρ, k])},PlotLabel-> "AutoCorrelation Function",PlotRange->All]

n