2.1 Modelos de média móvel (modelos MA) Modelos de séries temporais conhecidos como modelos ARIMA podem incluir termos auto-regressivos ou termos de média móvel. Na Semana 1, aprendemos um termo autorregressivo em um modelo de séries temporais para a variável x t é um valor retardado de x t. Por exemplo, um termo autorregressivo de atraso 1 é x t-1 (multiplicado por um coeficiente). Esta lição define termos de média móvel. Um termo de média móvel em um modelo de séries temporais é um erro passado (multiplicado por um coeficiente). Vamos (wt desviar N (0, sigma2w)), significando que os w t são identicamente, distribuídos independentemente, cada um com uma distribuição normal com média 0 e a mesma variância. O modelo de média móvel de ordem 1, denotado por MA (1) é (xt mu wt theta1w) O modelo de média móvel de 2ª ordem, denotado por MA (2) é (xt mu wt theta1w theta2w) , Denotado por MA (q) é (xt mu wt theta1w theta2w pontos thetaqw) Nota. Muitos livros didáticos e programas de software definem o modelo com sinais negativos antes dos termos. Isso não altera as propriedades teóricas gerais do modelo, embora ele inverta os sinais algébricos de valores de coeficientes estimados e de termos (não-quadrados) nas fórmulas para ACFs e variâncias. Você precisa verificar seu software para verificar se sinais negativos ou positivos foram usados para escrever corretamente o modelo estimado. R usa sinais positivos em seu modelo subjacente, como fazemos aqui. Propriedades Teóricas de uma Série de Tempo com um Modelo MA (1) Observe que o único valor não nulo na ACF teórica é para o atraso 1. Todas as outras autocorrelações são 0. Assim, uma ACF de amostra com uma autocorrelação significativa apenas no intervalo 1 é um indicador de um possível modelo MA (1). Para os estudantes interessados, provas destas propriedades são um apêndice a este folheto. Exemplo 1 Suponha que um modelo MA (1) seja x t 10 w t .7 w t-1. Onde (wt overset N (0,1)). Assim, o coeficiente 1 0,7. O ACF teórico é dado por Um gráfico deste ACF segue. O gráfico apenas mostrado é o ACF teórico para um MA (1) com 1 0,7. Na prática, uma amostra normalmente não proporciona um padrão tão claro. Usando R, simulamos n 100 valores de amostra usando o modelo x t 10 w t .7 w t-1 onde w t iid N (0,1). Para esta simulação, segue-se um gráfico de séries temporais dos dados da amostra. Não podemos dizer muito desse enredo. A ACF de amostra para os dados simulados segue. Observamos que a amostra ACF não corresponde ao padrão teórico do MA subjacente (1), ou seja, que todas as autocorrelações para os atrasos de 1 serão 0 Uma amostra diferente teria uma ACF de amostra ligeiramente diferente mostrada abaixo, mas provavelmente teria as mesmas características gerais. Propriedades teóricas de uma série temporal com um modelo MA (2) Para o modelo MA (2), as propriedades teóricas são as seguintes: Note que os únicos valores não nulos na ACF teórica são para os retornos 1 e 2. As autocorrelações para atrasos maiores são 0 . Assim, uma ACF de amostra com autocorrelações significativas nos intervalos 1 e 2, mas autocorrelações não significativas para atrasos maiores indica um possível modelo MA (2). Iid N (0,1). Os coeficientes são 1 0,5 e 2 0,3. Como este é um MA (2), o ACF teórico terá valores não nulos apenas nos intervalos 1 e 2. Os valores das duas autocorrelações não nulas são: Um gráfico do ACF teórico segue. Como quase sempre é o caso, dados de exemplo não vai se comportar tão perfeitamente como a teoria. Foram simulados n 150 valores de amostra para o modelo x t 10 w t .5 w t-1 .3 w t-2. Onde w t iid N (0,1). O gráfico de série de tempo dos dados segue. Como com o gráfico de série de tempo para os dados de amostra de MA (1), você não pode dizer muito dele. A ACF de amostra para os dados simulados segue. O padrão é típico para situações em que um modelo MA (2) pode ser útil. Existem dois picos estatisticamente significativos nos intervalos 1 e 2, seguidos por valores não significativos para outros desfasamentos. Note que devido ao erro de amostragem, a ACF da amostra não corresponde exactamente ao padrão teórico. ACF para Modelos Gerais MA (q) Uma propriedade dos modelos MA (q) em geral é que existem autocorrelações não nulas para os primeiros q lags e autocorrelações 0 para todos os retornos gt q. Não-unicidade de conexão entre os valores de 1 e (rho1) no modelo MA (1). No modelo MA (1), para qualquer valor de 1. O recíproco 1 1 dá o mesmo valor para Como exemplo, use 0,5 para 1. E então use 1 (0,5) 2 para 1. Você obterá (rho1) 0,4 em ambas as instâncias. Para satisfazer uma restrição teórica chamada invertibilidade. Restringimos modelos MA (1) para ter valores com valor absoluto menor que 1. No exemplo dado, 1 0,5 será um valor de parâmetro permitido, enquanto 1 10,5 2 não. Invertibilidade de modelos MA Um modelo MA é dito ser inversível se for algébrica equivalente a um modelo de ordem infinita convergente. Por convergência, queremos dizer que os coeficientes de RA diminuem para 0 à medida que avançamos no tempo. Invertibilidade é uma restrição programada em séries temporais de software utilizado para estimar os coeficientes de modelos com MA termos. Não é algo que verificamos na análise de dados. Informações adicionais sobre a restrição de invertibilidade para modelos MA (1) são fornecidas no apêndice. Teoria Avançada Nota. Para um modelo MA (q) com um ACF especificado, existe apenas um modelo invertible. A condição necessária para a invertibilidade é que os coeficientes têm valores tais que a equação 1- 1 y-. - q y q 0 tem soluções para y que caem fora do círculo unitário. Código R para os Exemplos No Exemplo 1, traçamos o ACF teórico do modelo x t 10w t. 7w t-1. E depois simularam n 150 valores a partir deste modelo e traçaram a amostra de séries temporais ea amostra ACF para os dados simulados. Os comandos R utilizados para traçar o ACF teórico foram: acfma1ARMAacf (mac (0.7), lag. max10) 10 lags de ACF para MA (1) com theta1 0.7 lags0: 10 cria uma variável chamada lags que varia de 0 a 10. plot (Lags, acfma1, xlimc (1,10), ylabr, typeh, ACF principal para MA (1) com theta1 0,7) abline (h0) adiciona um eixo horizontal ao gráfico O primeiro comando determina o ACF e o armazena em um objeto Chamado acfma1 (nossa escolha de nome). O comando de plotagem (o terceiro comando) traça defasagens em relação aos valores de ACF para os retornos de 1 a 10. O parâmetro ylab marca o eixo y eo parâmetro principal coloca um título no gráfico. Para ver os valores numéricos do ACF basta usar o comando acfma1. A simulação e as parcelas foram feitas com os seguintes comandos. Xcarima. sim (n150, lista (mac (0.7))) Simula n 150 valores de MA (1) xxc10 adiciona 10 para fazer a média 10. Padrões de simulação significam 0. plot (x, typeb, mainSimulated MA (1) data) Acf (x, xlimc (1,10), mainACF para dados de amostras simulados) No Exemplo 2, traçamos o ACF teórico do modelo xt 10 wt. 5 w t-1 .3 w t-2. E depois simularam n 150 valores a partir deste modelo e traçaram a amostra de séries temporais ea amostra ACF para os dados simulados. Os comandos R utilizados foram acfma2ARMAacf (mac (0,5,0,3), lag. max10) acfma2 lags0: 10 parcela (lags, acfma2, xlimc (1,10), ylabr, tipoh, ACF principal para MA (2) com theta1 0,5, (X, typeb, main Simulado MA (2) Series) acf (x, xlimc (1,10), x2, MainACF para dados simulados de MA (2) Apêndice: Prova de Propriedades de MA (1) Para estudantes interessados, aqui estão as provas para propriedades teóricas do modelo MA (1). Quando h 1, a expressão anterior 1 w 2. Para qualquer h 2, a expressão anterior 0 (x) é a expressão anterior x (x) A razão é que, por definição de independência do wt. E (w k w j) 0 para qualquer k j. Além disso, porque w t tem média 0, E (w j w j) E (w j 2) w 2. Para uma série de tempo, aplique este resultado para obter o ACF fornecido acima. Um modelo MA reversível é aquele que pode ser escrito como um modelo de ordem infinita AR que converge de modo que os coeficientes AR convergem para 0 à medida que nos movemos infinitamente para trás no tempo. Bem demonstrar invertibilidade para o modelo MA (1). Em seguida, substitui-se a relação (2) para wt-1 na equação (1) (3) (zt wt theta1 (z-theta1w) wt theta1z-theta2w) No tempo t-2. A equação (2) torna-se Então substituimos a relação (4) para wt-2 na equação (3) (zt wt theta1 z - theta21w wt theta1z - theta21 (z - theta1w) wt theta1z-theta12z theta31w) Se continuássemos Infinitamente), obteríamos o modelo AR de ordem infinita (zt wt theta1 z - theta21z theta31z - theta41z pontos) Observe, no entanto, que se 1 1, os coeficientes multiplicando os desfasamentos de z aumentarão (infinitamente) Tempo. Para evitar isso, precisamos de 1 lt1. Esta é a condição para um modelo MA (1) invertible. Infinite Order MA model Na semana 3, bem ver que um modelo AR (1) pode ser convertido em um modelo de ordem infinita MA: (xt - mu wt phi1w phi21w pontos phik1 w dots sum phij1w) Esta soma de termos de ruído branco passado é conhecido Como a representação causal de um AR (1). Em outras palavras, x t é um tipo especial de MA com um número infinito de termos voltando no tempo. Isso é chamado de ordem infinita MA ou MA (). Uma ordem finita MA é uma ordem infinita AR e qualquer ordem finita AR é uma ordem infinita MA. Lembre-se na Semana 1, observamos que um requisito para um AR estacionário (1) é que 1 lt1. Vamos calcular o Var (x t) usando a representação causal. Esta última etapa usa um fato básico sobre séries geométricas que requer (phi1lt1) caso contrário, a série diverge. NavigationDocumentation a é um vetor constante de deslocamentos, com n elementos. A i são n-by-n matrizes para cada i. Os Ai são matrizes autorregressivas. Existem p matrizes autorregressivas. 949 t é um vector de inovações não correlacionadas em série. Vetores de comprimento n. Os 949 t são vetores aleatórios normais multivariados com matriz de covariância Q. Onde Q é uma matriz de identidade, salvo indicação em contrário. B j são n-by-n matrizes para cada j. As B j são matrizes de média móvel. Existem q matrizes de média móvel. Xt é uma matriz n-by-r representando termos exógenos em cada momento t. R é o número de séries exógenas. Termos exógenos são dados (ou outras entradas não modificadas) além da série de tempo de resposta y t. B é um vetor constante de coeficientes de regressão de tamanho r. Portanto, o produto X t middotb é um vetor de tamanho n. Geralmente, as séries temporais y t e X t são observáveis. Em outras palavras, se você tiver dados, ele representa uma ou ambas as séries. Você nem sempre sabe o deslocamento a. Coeficiente b. Matrizes autorregressivas A i. E matrizes de média móvel B j. Normalmente, você deseja ajustar esses parâmetros aos seus dados. Consulte a página de referência da função vgxvarx para obter formas de estimar parâmetros desconhecidos. As inovações 949 t não são observáveis, pelo menos em dados, embora possam ser observadas em simulações. Lag Representação do Operador Há uma representação equivalente das equações auto-regressivas lineares em termos de operadores de retardamento. O operador de atraso L move o índice de tempo de volta em um: L y t y t 82111. O operador L m move o índice de tempo para trás por m. L m y t y t 8211 m. Na forma de operador de lag, a equação para um modelo de SVARMAX (p. Q. R) torna-se (A 0 x 2212 x2211 i 1 p A i L i) y t a X t b (B 0 x 2211 j 1 q B j L j) x03B5 t. Esta equação pode ser escrita como A (L) y t a X t b B (L) x03B5 t. Um modelo VAR é estável se det (I n x2212 A 1 z x 2212 A 2 z 2 x 2212. x2212 A pzp) x2260 0 x00A0x00A0forx00A0x00A0 z x2264 1. Esta condição implica que, com todas as inovações igual a zero, o processo VAR converge para um Como o tempo passa. Veja Luumltkepohl 74 Capítulo 2 para uma discussão. Um modelo VMA é inversível se det (I n B 1 z B 2 z 2. B q z q) x2260 0 x00A0x00A0forx00A0x00A0 z x2264 1. Esta condição implica que a representação VAR pura do processo é estável. Para obter uma explicação de como converter entre modelos VAR e VMA, consulte Alterando Representações de Modelo. Veja o capítulo 11 de Luumltkepohl para uma discussão sobre os modelos VMA invertíveis. Um modelo VARMA é estável se sua parte VAR é estável. Da mesma forma, um modelo VARMA é inversível se sua parte VMA é invertible. Não existe uma noção bem definida de estabilidade ou de invertibilidade para modelos com entradas exógenas (por exemplo, modelos VARMAX). Uma entrada exógena pode desestabilizar um modelo. Construindo Modelos VAR Para entender um modelo de séries temporais múltiplas, ou vários dados de séries temporais, geralmente você executa as seguintes etapas: Importe e pré-processa dados. Especifique um modelo. Estruturas de Especificação sem Valores de Parâmetro para especificar um modelo quando você deseja que o MATLAB x00AE estime os parâmetros Estruturas de Especificação com Valores de Parâmetro Selecionados para especificar um modelo onde você conhece alguns parâmetros e deseja que o MATLAB estime os outros Determinando um Número Apropriado de Lags a determinar Um número adequado de defasagens para seu modelo Ajustar o modelo aos dados. Ajustando Modelos a Dados para usar o vgxvarx para estimar os parâmetros desconhecidos em seus modelos. Isso pode envolver: Modificação de Representações de Modelo para alterar seu modelo para um tipo que vgxvarx manipula Analisar e prever usando o modelo ajustado. Isto pode envolver: Examinando a estabilidade de um modelo ajustado para determinar se seu modelo é estável e invertible. Modelo VAR Previsão para prever diretamente a partir de modelos ou para prever usando uma simulação de Monte Carlo. Cálculo de respostas de impulso para calcular as respostas de impulso, que fornecem previsões baseadas numa alteração assumida numa entrada para uma série temporal. Compare os resultados das previsões de seus modelos com os dados disponíveis para a previsão. Para um exemplo, veja Estudo de Caso do Modelo VAR. Seu aplicativo não precisa envolver todas as etapas deste fluxo de trabalho. Por exemplo, você pode não ter quaisquer dados, mas sim simular um modelo parametrizado. Nesse caso, você executaria apenas as etapas 2 e 4 do fluxo de trabalho genérico. Você pode iterar através de algumas dessas etapas. Exemplos Relacionados Selecione seu País11.2: Vector Modelos Autoregressivos Modelos VAR (p) Modelos VAR (modelos autorregressivos de vetor) são usados para séries temporais multivariadas. A estrutura é que cada variável é uma função linear de atrasos passados de si mesmo e atrasos passados das outras variáveis. Como exemplo, suponha que medimos três variáveis de séries temporais diferentes, denotadas por (x), (x) e (x). O modelo vectorial autorregressivo de ordem 1, denotado como VAR (1), é o seguinte: Cada variável é uma função linear dos valores de atraso 1 para todas as variáveis no conjunto. Em um modelo VAR (2), os valores de atraso 2 para todas as variáveis são adicionados ao lado direito das equações. No caso de três variáveis x (ou séries temporais) haveria seis preditores no lado direito de cada equação , Três lag 1 termos e três lag 2 termos. Em geral, para um modelo VAR (p), os primeiros p lags de cada variável no sistema seriam usados como preditores de regressão para cada variável. Os modelos VAR são um caso específico de modelos VARMA mais gerais. Os modelos VARMA para séries temporais multivariadas incluem a estrutura VAR acima, juntamente com os termos de média móvel para cada variável. Mais geralmente ainda, estes são casos especiais de modelos ARMAX que permitem a adição de outros preditores que estão fora do conjunto multivariado de interesse principal. Aqui, como na Seção 5.8 do texto, bem focar modelos VAR. Na página 304, os autores se encaixam no modelo da forma mathbf t Gamma mathbf t phi matemática mathbf t onde (mathbf t (1, t)) inclui termos para ajustar simultaneamente a constante ea tendência. Ele surgiu a partir de dados macroeconômicos onde grandes mudanças nos dados permanentemente afetam o nível da série. Há uma diferença não tão sutil aqui das lições anteriores em que agora estamos ajustando um modelo a dados que não precisam ser estacionários. Em versões anteriores do texto, os autores separadamente de-tendência de cada série usando uma regressão linear com t, o índice de tempo, como a variável preditor. Os valores de desvios para cada uma das três séries são os resíduos dessa regressão linear em t. O desvirtuamento é útil conceitualmente porque tira a força de direção comum que o tempo pode ter em cada série e criou stationarity como nós vimos em lições passadas. Esta abordagem resulta em coeficientes semelhantes, embora ligeiramente diferentes, uma vez que agora estamos simultaneamente ajustando a intercepção ea tendência juntos em um modelo OLS multivariável. A biblioteca de Vars authored by Bernhard Pfaff tem a capacidade de se encaixar neste modelo com tendência. Vejamos dois exemplos: um modelo estacionário diferencial e um modelo estacionário-tendencial. Modelo Diferencial-Estacionário O Exemplo 5.10 do texto é um modelo diferencial-estacionário em que as primeiras diferenças são estacionárias. Se examinarmos o código e o exemplo do texto ajustando o modelo acima: install. packages (vars) Se ainda não estiver instalado install. packages (astsa) Se não estiver instalado biblioteca (vars) biblioteca (astsa) x cbind (cmort, Parte) plot. ts (x, main, xlab) sumário (VAR (x, p1, typeboth)) Os dois primeiros comandos carregam os comandos necessários da biblioteca vars e os dados necessários da nossa biblioteca de textos. O comando cbind cria um vetor de variáveis de resposta (um passo necessário para respostas multivariadas). O comando VAR faz estimativa de modelos AR usando mínimos quadrados ordinários enquanto simultaneamente ajusta a tendência, o intercepto e o modelo ARIMA. O argumento p 1 solicita uma estrutura AR (1) e ambos se encaixam constante e tendência. Com o vetor de respostas, seu realmente um VAR (1). A seguir está a saída do comando VAR para a variável tempr (o texto fornece a saída para cmort): Os coeficientes para uma variável são listados na coluna Estimativa. O. l1 anexado a cada nome de variável indica que elas são variáveis de atraso 1. Usando a notação T temperatura, ttime (coletado semanalmente), taxa de mortalidade M, e P de poluição, a equação para a temperatura é de 67.586 - .007 t - 0.244 M 0.487 T - 0.128 P A equação para a taxa de mortalidade é de 73.227 t 0.014 t 0,465 M - 0,361 T 0,099 P A equação para a poluição é de 67,464 - 0,005 t - 0,125 M - 0,477 T 0,581 P. A matriz de covariância dos resíduos do VAR (1) para as três variáveis é impressa abaixo dos resultados da estimação. As variâncias estão abaixo da diagonal e podem ser usadas para comparar este modelo com VARs de ordem superior. O determinante dessa matriz é usado no cálculo da estatística BIC que pode ser usada para comparar o ajuste do modelo ao ajuste de outros modelos (ver fórmulas 5.89 e 5.90 do texto). Para mais referências sobre esta técnica ver Análise de séries temporais integradas e co-integradas com R por Pfaff e também Campbell e Perron 1991. No Exemplo 5.11 na página 307, os autores apresentam resultados para um modelo VAR (2) para os dados de taxa de mortalidade . Em R, você pode ajustar o modelo VAR (2) com o resumo do comando (VAR (x, p2, typeboth)) A saída, conforme exibido pelo comando VAR é a seguinte: Novamente, os coeficientes para uma variável particular são listados em A coluna Estimativa. Como exemplo, a equação estimada para a temperatura é a seguinte: 49.88 - .005 t - 0.109 M 0.261 T 0.051 P - 0.041 M 0.356 T 0.095 P Discutiremos estatísticas de critério de informação para comparar modelos VAR de diferentes ordens no trabalho de casa. Os resíduos também estão disponíveis para análise. Por exemplo, se atribuirmos o comando VAR a um objeto intitulado fitvar2 no nosso programa, fitvar2 VAR (x, p2, typeboth) então temos acesso aos resíduos da matriz (fitvar2). Esta matriz terá três colunas, uma coluna de resíduos para cada variável. Por exemplo, podemos usar para ver a ACF dos resíduos para a taxa de mortalidade depois de ajustar o modelo VAR (2). A seguir está o ACF que resultou do comando que acabamos de descrever. Parece bom para um ACF residual. (O pico grande no início é a correlação de atraso 0 sem importância.) Os dois comandos a seguir criarão ACFs para os resíduos para as outras duas variáveis. Eles também se assemelham a ruído branco. Podemos também examinar esses gráficos na matriz de correlação cruzada fornecida por acf (resíduos (fitvar2)): As parcelas ao longo da diagonal são as ACFs individuais para cada modelo de resíduos que acabamos de discutir acima. Além disso, vemos agora os gráficos de correlação cruzada de cada conjunto de resíduos. Idealmente, estes também se assemelham a ruído branco, no entanto, vemos restante cross-correlações, especialmente entre a temperatura ea poluição. Como observam nossos autores, este modelo não capta adequadamente a associação completa entre essas variáveis no tempo. Modelo Trend-Estacionário Permite explorar um exemplo onde os dados originais estão estacionários e examinar o código VAR, ajustando o modelo acima com uma constante e tendência. Usando R, simulamos n 500 valores de amostra usando o modelo VAR (2) Usando o comando VAR explicado acima: y1scan (var2daty1.dat) y2scan (var2daty2.dat) resumo (VAR (cbind (y1, y2), p2, typeboth) ) Obtemos o seguinte resultado: As estimativas são muito próximas dos coeficientes simulados ea tendência não é significativa, como esperado. Para dados estacionários, quando a destruição é desnecessária, você também pode usar o comando ar. ols para ajustar um modelo VAR: fitvar2 ar. ols (cbind (y1, y2), order2) Na primeira matriz dada, ler em uma linha para obter Os coeficientes para uma variável. As vírgulas precedentes, seguidas de 1 ou 2, indicam se os coeficientes são variáveis de atraso 1 ou atraso 2, respectivamente. As intercepções das equações são dadas em x. intercept um intercepto por variável. A matriz sob var. pred dá a matriz de variância-covariância dos resíduos da VAR (2) para as duas variáveis. As variâncias estão abaixo da diagonal e podem ser usadas para comparar este modelo com VARs de ordem superior como observado acima. Os erros padrão dos coeficientes AR são dados pelo comando fitvar2asy. se. coef. A saída é Como com os coeficientes, leia em linhas. A primeira linha dá os erros-padrão dos coeficientes para as variáveis de atraso 1 que predizem y1. A segunda linha dá os erros-padrão para os coeficientes que predizem y2. Você pode notar que os coeficientes estão próximos ao comando VAR, exceto o intercepto. Isso ocorre porque ar. ols estima o modelo para x-mean (x). Para coincidir com a intercepção fornecida pelo comando resumo (VAR (cbind (y1, y2), p2, typeconst)), você deve calcular a intercepção da seguinte maneira: Em nosso exemplo, a interceptação para o modelo simulado para yt, 1 é igual a -0,043637 -2.733607 (1-0.29300.4523) 15.45479 (-0.1913-0.6365) 9.580768 e a equação estimada para yt, 1 Estimativa com Minitab Para usuários do Minitab, está o fluxo geral do que fazer. Leia os dados em colunas. Utilize as séries temporais gt Lag para criar as colunas atrasadas necessárias dos valores estacionários. Utilizar o Stat gt ANOVA gt Geral MANOVA. Insira a lista de variáveis de tempo presentes como variáveis de resposta. Insira as variáveis x defasadas como covariáveis (e como o modelo). Clique em Resultados e selecione Análise Univariada (para ver os coeficientes de regressão estimados para cada equação). Se desejar, clique em Armazenamento e selecione Resíduos e / ou Ajustes. NavigationMoving-Average Representação de Aproximações Autoregressivas Nós estudamos as propriedades de uma infinita MA-representação de uma aproximação autorregressiva para um estacionário, real-valorizado processo. Ao fazer isso, damos uma extensão do Teorema de Wieners na configuração de aproximação determinística. Ao lidar com dados, podemos usar este novo resultado chave para obter insights sobre a estrutura de infinitas MA-representações de modelos auto-regressivos ajustados onde a ordem aumenta com o tamanho da amostra. Em particular, nós damos um limite uniforme para estimar os coeficientes de média móvel através de aproximação autorregressiva uniforme em todos os inteiros. 423.pdf
No comments:
Post a Comment